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I.  INTRODUCTION 

The  survivability  of  naval  ships  and  submarines  is  of  supreme  importance  to  those  who 
design,  build  and  sail  them.  Underwater  explosions,  created  by  the  detonation  of  mines  or 
torpedoes  near  naval  vessels,  clearly  represent  a  significant  threat  to  that  survivability.  The 
capacity  of  such  weapons  to  cause  major  damage  to  naval  vessels  has  been  recognized  since 
the  beginning  of  the  age  of  modern  naval  warfare,  and  has  been  the  subject  of  extensive 
investigation  since  the  earliest  years  of  the  Second  World  War. 

One  result  of  the  allied  experience  during  the  Second  World  War  was  the  realization 
that  the  pressure  wave  resulting  from  the  pulsation  (oscillation)  of  the  bubble  of  gases 
produced  by  the  detonation  of  an  explosive  material  underwater  can  be  as  important  as  the 
primary  shock  wave  in  causing  damage  to  a  naval  vessel.  As  Hicks  (1970a)  has  noted,  this 
bubble  pulsation  is  the  primary  cause  of  the  observed  "whipping"  of  ships,  and  often  has  a 
period  near  that  of  the  hull  girder's  fundamental  flexural  vibration  mode. 

This  will,  in  the  worst  case,  break  the  back  of  the  ship.  Even  in  less  severe  situations, 
the  oscillation  of  the  ship's  hull  may  be  significantly  increased  by  the  timing  of  the  arrival  of 
the  secondary  pressure  pulse  radiated  when  the  bubble  collapses  to  its  first  minimum  volume. 
This  can  cause  appreciable  damage  to  sensitive  internal  equipment.  Furthermore,  the 
increased  use  of  commercial  off-the-shelf  (COTS)  equipment  within  modern  ships,  while 
fiscally  necessary  and  undoubtedly  beneficial  in  the  long  run,  may  increase  the  damage 
potential  of  explosion  gas  bubble  oscillations. 

Fundamental  analysis  shows  that  low  frequency  vibration  input  can  significantly 
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increase  the  resulting  response  in  certain  underdamped  system  configurations.  Equipment  not 
designed  to  withstand  the  increased  response  could  easily  suffer  damage  that,  depending  upon 
ship  design,  reduces  the  war  fighting  capability  of  the  ship.  In  an  age  of  over-the-horizon 
saturation  targeting,  such  a  calamity  may  be  as  bad  as  breaking  the  back  of  the  ship. 

A  significant  investment  of  manpower  and  time  has  been  expended  in  analyzing  the 
phenomena  associated  with  the  pulsation  of  explosion  gas  bubbles,  both  during  and  since  the 
Second  World  War.  However,  these  studies  have,  for  the  most  part,  relied  upon  many 
simplifying  assumptions.  The  validity  of  some  of  these  assumptions  were  known  beforehand 
to  be  very  inaccurate  during  certain  phases  of  an  explosion  gas  bubble's  motion.  In  fact,  some 
of  these  assumptions  are  least  accurate  during  the  times  when  key  explosion  gas  bubble 
phenomena  are  occurring.  Predictions  of  explosion  gas  bubble  behavior  from  these  studies 
are  therefore  more  likely  to  be  qualitatively  correct  than  quantitatively  accurate.  This  is  an 
important  consideration  because  ship  design  decisions  are  based,  in  part,  upon  these 
predictions. 

The  next  section  describes  some  of  the  important  phenomena  associated  with 
underwater  explosion  gas  bubbles.  Following  this  is  a  synopsis  of  the  more  significant 
research  that  has  been  previously  conducted  into  underwater  explosion  gas  bubbles.  The  third 
and  final  section  of  this  introductory  chapter  discusses  the  objectives  of  this  present  research 
effort,  and  presents  an  overview  of  the  remaining  chapters. 
A.         BUBBLE  PHENOMENA 

Consider  a  spherical  explosive  charge  located  underwater  at  a  considerable  distance 
from  any  boundary  surface.  When  this  charge  is  detonated,  the  explosive  material  is  very 


rapidly  converted  into  high  pressure  gaseous  reaction  products,  and  a  high  pressure  shock 
wave  is  transmitted  to  and  propagated  through  the  surrounding  fluid.  About  50%  of  the 
initial  chemical  energy  of  the  explosive  is  transmitted  to  the  fluid  in  this  initial  shock  wave 
(Snay,  1957). 

Because  this  shock  wave  propagates  radially  outward,  the  amplitude  of  the  shock 
wave  decreases  with  increasing  distance  from  the  center  of  the  charge.  The  fluid  behind  the 
shock  wave  front  attains  a  large  outward  radial  velocity.  This  causes  the  pressure  at  some 
distance  behind  the  shock  wave  front  to  drop  below  the  hydrostatic  pressure  of  the  fluid 
ahead  of  this  wave  front. 

The  pressure  in  the  gas  bubble  is  significantly  reduced  after  emission  of  the  primary 
shock  wave,  but  is  still  significantly  higher  than  the  hydrostatic  pressure  in  the  surrounding 
fluid.  This,  in  conjunction  with  the  radial  outflow  of  fluid  from  the  vicinity  of  the  bubble, 
causes  the  bubble  to  expand  rapidly.  This  expansion  continues  for  a  relatively  long  time,  the 
pressure  in  the  bubble  decreasing  as  the  bubble  volume  increases.  The  pressure  in  the  gas 
bubble  eventually  falls  below  hydrostatic  pressure,  but  the  expansion  persists  because  of  the 
inertia  of  the  outward  flowing  fluid  (Cole,  1948,  p.  8). 

Eventually,  the  fluid  immediately  surrounding  the  bubble  comes  to  rest.  This,  together 
with  the  very  low  pressure  in  the  bubble  at  this  time,  causes  the  bubble  to  begin  contracting. 
The  fluid  in  the  immediate  vicinity  of  the  bubble  begins  flowing  inward,  accelerating  the 
contraction  of  the  bubble.  As  the  volume  of  the  bubble  decreases,  its  internal  pressure 
increases.  The  pressure  within  the  gas  bubble  eventually  becomes  significantly  higher  than 
the  hydrostatic  pressure  of  the  surrounding  fluid,  and  the  contraction  of  the  bubble  is  reversed 


abruptly.  The  inertia  of  the  surrounding  fluid  together  with  the  compressibility  of  the  bubble 
gases  and  the  surrounding  fluid  "thus  provide  the  necessary  conditions  for  an  oscillating 
system,  and  the  bubble  does  in  fact  undergo  repeated  cycles  of  expansion  and  contraction" 
(Cole,  1948,  p.  8). 

Figure  1-1  illustrates  the  bubble  radius  vs.  time  behavior  typically  observed  in  an 
underwater  explosion.  Note  that  the  maximum  radius  of  the  bubble  decreases  during 
subsequent  expansions.  This  is  due,  at  least  in  part,  to  the  emission  of  secondary  pressure 
pulses  when  the  bubble  radius  is  near  its  minimum  values.  About  66%  of  the  energy 
remaining  in  the  gas  bubble  after  emission  of  the  primary  shock  wave  is  lost  during  the  first 
expansion-contraction  cycle  (Cole,  1948,  p.  283).  Successive  bubble  pulses  are  thus  weaker, 
and  generally  only  the  first  pulse  is  of  practical  significance  (Cole,  1948,  p.  10). 

Figure  1-2  shows  a  typical  pressure  vs.  time  curve  at  fixed  distance  from  an 
underwater  explosion.  Although  the  peak  pressure  in  the  secondary  pressure  pulse  is  much 
less  than  the  peak  shock  wave  pressure,  the  duration  of  this  pulse  is  greater.  The  impulses 
due  to  the  primary  shock  wave  and  the  first  of  the  bubble  pulses  are  typically  comparable 
(Cole,  1948,  p.  364). 

Additional  phenomena  are  observed  when  an  explosion  occurs  near  a  boundary 
surface,  such  as  the  surface  or  bottom  of  the  ocean  or  a  nearby  structure.  When  an  explosion 
occurs  near  a  relatively  rigid  boundary  such  as  the  hull  of  a  naval  vessel  or  a  hard  ocean 
bottom,  the  explosion  gas  bubble  migrates  toward  the  boundary  during  its  contraction.  The 
opposite  effect  is  seen  when  an  explosion  occurs  near  a  free  surface  (Snay,  1957).  The  shape 
of  migrating  gas  bubbles  has  also  been  seen  to  be  non-spherical,  as  illustrated  in  Figure  1-3. 


Figure  1-1.  Typical  Bubble  Radius  vs.  Time  Curve  [From  Swift  and  Decius  (1950)] 
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Figure  1-2.      Typical  Pressure  Profile  Produced  From  an  Underwater  Explosion  [From 
Cole  (1948)] 


Figure  1-3.      Typical  Shape  for  a  Migrating  Underwater  Explosion  Gas  Bubble  [From 
Cole  (1948)] 


Since  the  secondary  pressure  pulse  is  emitted  when  a  bubble  is  near  its  minimum 
volume,  the  migration  of  a  bubble  determines  where  this  bubble  pulse  is  emitted  from.  For 
an  explosion  occurring  near  a  naval  ship,  the  migration  of  the  bubble  can  cause  the  bubble 
pulse  to  occur  much  closer  to  the  hull,  potentially  increasing  the  damage  caused  to  the  ship 
(Hicks,  1970b). 

The  presence  of  boundary  surfaces  also  affects  the  period  of  the  oscillation  of  the 
bubble.  A  rigid  boundary  causes  the  period  of  oscillation  of  the  bubble  to  increase,  and  a  free 
surface  causes  it  to  decrease  (Herring,  1950).  Thus,  if  an  explosion  occurs  near  the  hull  of 
a  ship,  the  presence  of  the  ship  itself  can  cause  the  time  interval  between  arrival  of  the  primary 
shock  wave  and  secondary  pressure  pulse  to  become  either  closer  to  or  further  from  the 
fundamental  period  of  the  hull  girder,  significantly  affecting  the  amount  of  whipping 
undergone  by  the  ship. 
B.         PREVIOUS  RESEARCH  INTO  EXPLOSION  GAS  BUBBLES 

Perhaps  the  earliest  contribution  related  directly  to  an  appreciation  and  understanding 
of  bubbles  was  made  by  Reynolds  (1894),  who  noted  the  formation  of  vapor  cavitation 
bubbles  in  water  flowing  through  a  constricted  pipe.  The  first  analysis  of  the  dynamic 
behavior  of  bubbles  was  made  by  Rayleigh  (1917),  who  used  conservation  of  momentum  to 
derive  an  equation  for  the  collapse  of  a  spherical  void. 

Lamb  (1923)  studied  the  expansion  phase  of  a  gas  bubble,  and  derived  an  expression 
relating  the  maximum  radius  of  a  bubble  to  the  depth  and  total  energy,  neglecting  fluid 
compressibility.  Ramsauer  (1923)  conducted  small  scale  experiments  with  guncotton  which 
showed  reasonably  good  agreement  with  Lamb's  predicted  relationship  between  maximum 
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bubble  radius  and  depth.  The  work  of  Lamb  is  the  basis  for  the  standard  semi-empirical 
scaling  formula  for  the  maximum  radius  of  an  explosion  gas  bubble  when  buoyancy  can  be 
neglected  and  no  boundary  surfaces  are  near  the  charge.  This  situation  is  referred  to  as  the 
"free-field"  case.  This  scaling  formula  defines  a  relationship  between  the  maximum  radius  of 
the  bubble,  the  charge  weight,  and  the  charge  depth.  It  is  termed  "semi-empirical"  because 
it  depends  upon  a  proportionality  constant  which  varies  between  different  types  of  explosives, 
and  must  be  determined  through  experimentation. 

Willis  (1941)  was  the  first  to  succeed  in  integrating  the  energy  equation  for  the 
noncompressive  radial  motion  of  free-field  explosion  gas  bubbles,  and  thereby  derive  an 
expression  for  the  period  of  the  first  oscillation  cycle  of  the  bubble.  His  work  is  the  basis  for 
the  standard  semi-empirical  formula  for  the  free-field  first  bubble  oscillation  period  of  an 
explosion  gas  bubble.  This  formula  defines  a  scaling  relationship  for  different  charge  weights 
and  depths,  and  depends  upon  an  experimentally  determined  proportionality  constant.  Willis 
also  conducted  experiments  to  confirm  the  predicted  relationship  between  charge  size,  depth, 
and  bubble  period,  and  found  good  agreement.  This  result  has  since  been  confirmed  by  other 
researchers,  particularly  as  to  the  functional  dependency  of  the  first  bubble  period  on  the 
charge  depth  (Cole,  1948,  pp.  280-281). 

Herring  (1941,  1950)  made  a  number  of  important  contributions  to  the  current 
understanding  of  explosion  gas  bubble  phenomena.  He  predicted,  based  upon  an  approximate 
theoretical  development,  that  one  of  the  effects  of  fluid  compressibility  would  be  to  make  the 
bubble  radius  versus  time  curve  for  a  free-field  bubble  asymmetrical,  with  the  contraction 
being  slower  the  expansion.  He  conducted  an  analysis  indicating  that  the  radiation  of  energy, 


if  fluid  compressibility  was  allowed  for,  would  be  negligible  when  the  bubble  is  large,  but  that 
an  appreciable  radiation  of  energy  could  take  place  at  the  first  bubble  minimum. 

Herring's  most  important  contribution  was  the  advancement  of  an  approximate  theory 
to  account  for  the  effect  of  simple  rigid  boundaries  or  free  surfaces  (constant  pressure 
boundaries)  on  explosion  gas  bubbles  during  the  first  oscillation  cycle.  In  this  theory,  the 
effects  of  the  boundary  surface  are  treated  as  small  perturbations  of  the  motion  of  the  gas 
bubble  in  the  absence  of  the  boundary.  The  boundary  is  assumed  to  be  remote  enough  so  that 
the  standoff  distance  to  the  boundary  is  large  compared  to  the  maximum  radius  of  the  bubble, 
which  is  assumed  to  remain  nearly  spherical.  This  theory  is  thus  only  potentially  valid  when 
the  boundary  surface  is  fairly  remote  from  the  gas  bubble.  Furthermore,  the  theory  neglects 
fluid  compressibility  entirely. 

One  of  the  effects  of  fluid  compressibility  is  a  finite  wave  speed  in  the  fluid;  in  an 
incompressible  fluid,  the  wave  velocity  is  infinite.  Thus  in  a  compressible  fluid  there  is  a  finite 
distance  beyond  which  the  presence  of  a  boundary  cannot  affect  an  explosion  gas  bubble 
during  the  first  oscillation  cycle,  because  even  the  primary  shock  wave  resulting  from 
detonation  of  the  charge  will  not  have  reached  the  boundary  by  the  end  of  the  first  cycle.  The 
theory  developed  by  Herring  to  account  for  simple  boundary  surfaces  therefore  cannot  be 
correct  at  large  standoff  distances,  as  it  assumes  an  incompressible  fluid. 

Although  Herring's  theory  for  predicting  the  effect  of  these  simple  types  of  boundary 
surfaces  is  still  in  use,  e.g.  in  the  computer  code  MSWHlP  (Hicks,  1971)  (which  is  still  being 
used  by  the  United  States  Navy  to  predict  the  response  of  naval  vessels  to  underwater 
explosion  gas  bubbles),  the  above  discussion  begs  an  answer  to  the  question  of  when  it  is 
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valid.  It  was  never  intended  to  be  accurate  at  small  standoff  distances,  and  because  it  neglects 
fluid  compressibility,  it  is  not  valid  at  large  standoff  distances. 

Another  important  contribution  to  the  current  understanding  of  explosion  gas  bubble 
phenomena  was  made  by  Taylor  (1943),  who  derived  an  expression  for  the  vertical  migration 
of  a  spherically  symmetric  bubble  due  to  gravity  (buoyancy)  in  the  absence  of  nearby 
boundary  surfaces,  neglecting  fluid  compressibility.  In  addition  to  neglecting  fluid 
compressibility,  Taylor's  theory  is  based  upon  the  assumption  that  artificial  forces,  which  do 
no  work,  act  to  keep  the  gas  bubble  spherical.  Actual  experimental  results  typically  show 
significant  bubble  shape  departure  from  spherical,  as  illustrated  in  Figure  1-3.  In  the  present 
application  of  Taylor's  theory  for  the  effects  of  buoyancy  on  explosion  gas  bubbles,  use  is  also 
made  of  Herring's  theory  for  the  effect  of  the  free  surface  (if  a  free  surface  is  not  nearby,  the 
bubble  is  generally  deep  enough  so  that  buoyancy  is  not  significant). 

It  is  perhaps  then  not  surprising  that  predictions  of  explosion  gas  bubble  behavior  from 
the  combined  Herring-Taylor  theory  do  not  agree  very  well  with  experimental  results.  For 
example,  an  experimental  measurement  by  Bryant  (1950)  near  a  free  surface  showed  a  peak 
bubble  migration  velocity  of  180  ft/s  (55  m/s),  as  compared  to  a  calculated  value  of  890  ft/s 
(271  m/s). 

Numerous  researchers  have  conducted  underwater  explosion  experiments.  Notable 
among  these  is  Swift  and  Decius  (1950),  who  conducted  a  series  of  experiments  to  determine 
relatively  accurate  values  for  the  proportionality  constants  in  the  semi-empirical  formulas  for 
the  free-field  maximum  bubble  radius  and  first  bubble  period,  for  various  types  of  explosives. 

Chertock  (1952)  developed  an  approximate  method  to  calculate  the  bubble  pulse 
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induced  whipping  response  of  naval  ships  and  submarines,  for  explosions  occurring  very 
remote  from  the  hull.  Hicks  (1970a,  1970b,  1971)  determined  that  Chertock's  methodology 
might  be  applied  at  smaller  standoff  distances  with  simple  modifications,  and  then  further 
modified  this  procedure  to  account  for  bubble  migration,  using  the  Herring-Taylor  theory. 
In  an  attempt  to  account  for  the  known  problems  with  this  theory,  his  procedure  introduces 
an  artificially  large  drag  coefficient,  the  value  of  which  he  has  chosen  based  upon  very  limited 
experimental  data.  This  is  justified  as  being  "the  best  [solution]  available  until  fresh 
experimental  results  (particularly  for  the  near  field)  or  a  more  rational  theory  become 
available"  (Hicks,  1970b).  The  computer  code  developed  by  Hicks,  MSWHTP,  thus  far 
appears  to  generally  give  reasonably  accurate  results  for  distant  charges;  the  validity  of  its 
predictions  at  short  and  intermediate  standoff  distances  is  still  largely  unknown. 

Some  recent  efforts  have  examined  methods  to  study  underwater  explosion  gas  bubble 
phenomena  under  laboratory  conditions.  Schmidt  et.  al  (1987)  have  studied  the  problem  of 
conducting  very  small  scale  explosive  testing  (charge  weights  of  about  0.2  g)  while 
maintaining  complete  similarity,  by  using  a  centrifuge  to  obtain  gravitational  accelerations  of 
up  to  500  times  the  standard  terrestrial  gravitational  acceleration.  Chahine  et.  al.  (1995)  have 
developed  a  procedure  for  generating  plasma  bubbles  using  a  very  high  energy  spark  and 
scaling  these  plasma  bubbles  to  represent  explosion  gas  bubbles. 

As  can  be  seen  from  the  brief  summary  of  previous  work  given  above,  much  of  the 
current  understanding  and  theoretical  development  regarding  underwater  explosions  was  a 
result  of  allied  research  efforts  during  the  Second  World  War.  Many  of  the  more  important 
papers  from  this  era  can  be  found  in  the  three  volume  work  Underwater  Explosion  Research: 
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A  Compendium  of  British  and  American  Reports  (Office  of  Naval  Research,  1950).    A 
comprehensive  summary  of  the  research  from  this  period  is  provided  by  Cole  (1948). 
C.         OBJECTIVES  OF  CURRENT  RESEARCH 

The  above  synopsis  of  previous  research  into  explosion  gas  bubble  phenomena  points 
out  several  areas  in  which  the  existing  understanding  of  these  phenomena  is  inadequate.  The 
existing  theoretical  understanding  of  the  dynamic  behavior  of  explosion  gas  bubbles  is  based 
upon  neglecting  fluid  compressibility  entirely.  Only  limited,  qualitative  estimates  of  what  the 
affects  of  fluid  compressibility  might  be  have  been  made.  And  in  addition  to  neglecting  fluid 
compressibility,  the  existing  theory  describing  the  dynamic  behavior  of  explosion  gas  bubbles 
near  simple  rigid  or  constant  pressure  boundary  surfaces  does  not  cover  the  situation  in  which 
the  bubble  is  fairly  close  to  the  boundary.  Even  for  the  case  in  which  the  bubble  is  not  too 
close  to  the  boundary  surface,  the  dynamic  behavior  of  the  bubble  predicted  using  the  existing 
theory  is  known  to  be  more  qualitatively  correct  than  quantitatively  accurate. 

The  research  described  in  this  study  was  undertaken  to  rectify  some  of  these 
shortcomings  in  the  current  understanding  of  explosion  gas  bubble  phenomena.  It  was 
believed  that,  by  using  a  modern  numerical  analysis  computer  program  based  upon  a  finite 
control  volume  (Eulerian)  method,  it  would  be  possible  to  quantitatively  investigate  explosion 
gas  bubble  phenomena  which  had  heretofore  been  understood  only  qualitatively. 

The  first  goal  of  this  study  was  to  examine  the  feasibility  of  using  a  finite  control 
volume  numerical  analysis  technique  to  investigate  underwater  explosion  phenomena,  and  if 
this  method  proved  serviceable,  to  investigate  the  effects  of  fluid  compressibility  and  internal 
gas  energy  on  the  dynamic  behavior  of  explosion  gas  bubbles.  The  second  objective  was  to 
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investigate  the  effects  of  plane  rigid  and  constant  pressure  boundaries  on  important  aspects 
of  the  dynamic  behavior  of  explosion  gas  bubbles,  including  fluid  compressibility  and  internal 
energy,  and  to  quantitatively  characterize  these  effects. 

Two  specific  aspects  of  the  dynamic  behavior  of  bubbles  near  boundaries  were 
investigated.  The  first  of  these  was  the  change  in  the  period  of  the  bubble  caused  by  the 
boundary.  The  bubble  period  is  important  because  the  time  interval  between  arrival  of  the 
primary  shock  wave  and  the  first  bubble  pulse  can  be  very  near  the  fundamental  bending  mode 
period  of  a  ship's  hull  girder.  Thus,  if  the  boundary  causes  this  time  interval  to  become  even 
closer  to  this  hull  girder  bending  mode  period,  the  resultant  whipping  can  be  increased  due 
to  resonance.  As  discussed  earlier,  the  impulses  due  to  the  primary  shock  wave  and  the  first 
bubble  pulse  are  comparable  in  magnitude,  while  the  impulses  from  subsequent  bubble  pulses 
are  much  smaller.  The  other  aspect  of  a  bubble's  dynamic  behavior  that  was  investigated  was 
the  migration  of  the  bubble  at  the  end  of  the  first  expansion-contraction  cycle,  either  towards 
or  away  from  the  boundary.  Because  pressures  fall  off  with  increasing  distance,  it  is 
important  to  know  where  the  bubble  is  located  when  the  first  bubble  pulse  is  emitted. 

The  numerical  analysis  approach  used  for  this  research  is  described  in  the  next  chapter. 
Chapter  HI  describes  the  application  of  this  analysis  approach  to  an  investigation  of  the  affects 
of  fluid  compressibility  and  internal  gas  energy  on  a  free-field  explosion  gas  bubble,  and 
compares  results  from  this  analysis  approach  with  both  experimental  and  analytical  results. 
The  effects  of  rigid  and  constant  pressure  boundaries  on  explosion  gas  bubbles,  including  fluid 
compressibility  and  internal  gas  energy,  are  examined  in  Chapter  IV.  Chapter  V  then 
summarizes  the  results  from  this  study. 
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This  dissertation  illustrates  the  application  of  a  finite  volume  based  numerical  analysis 
technique  for  investigation  of  underwater  explosion  gas  bubbles.  This  technique  was  used  to 
investigate  the  effects  of  fluid  compressibility  and  internal  energy  on  a  free-field  explosion  gas 
bubble.  This  analysis  procedure  was  also  used  to  quantitatively  characterize  the  effects  of 
rigid  and  constant  pressure  boundaries  on  explosion  gas  bubbles,  when  internal  energy  and 
fluid  compressibility  are  not  neglected.  The  quantitative  characterizations  derived  have  a 
form  suitable  for  use  in  scaling  these  results  for  other  explosive  charge  types,  weights,  and 
depths. 
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H.  NUMERICAL  SOLUTION  TECHNIQUE 

Because  underwater  explosions  involve  the  flow  of  liquids  and  gases,  and  a  high 
pressure  shock  wave,  numerical  analysis  of  underwater  explosion  phenomena  in  the  region 
near  the  charge  using  standard  Lagrangian  based  finite  element  programs  is  not  practical.  In 
a  Lagrangian  based  finite  element  program,  the  elements  deform  in  response  to  the  pressure 
on  the  faces  of  the  element.  This  leads  to  elements  in  the  region  of  the  shock  front  being 
severely  deformed  (crushed),  and  the  time  step  size  per  iteration  becoming  very  small,  so  that 
the  solution  advances  very  slowly  in  time. 

The  numerical  analyses  described  in  this  study  were  therefore  conducted  using  an 
Eulerian  based  finite  volume  program.  The  computer  program  MSC/DYTRAN  (The 
MacNeal  Schwendler  Corporation,  1995)  was  used  for  the  numerical  analyses  conducted 
during  this  investigation.  The  analyses  used  the  multi-material  Eulerian  processor  in  this 
program,  which  is  based  upon  the  computer  program  MSC/PISCES  (The  MacNeal 
Schwendler  Corporation,  1991).  This  processor  provides  the  flexibility  of  general 
connectivity  for  Eulerian  finite  volumes,  so  the  finite  volumes  are  not  restricted  as  to  shape. 
This  feature  is  useful  when  a  large  volume  of  fluid  must  be  discretized,  but  the  area  of  interest 
is  relatively  small.  Other  Eulerian  based  programs  typically  require  a  uniform  discretization, 
which,  for  analysis  of  an  underwater  explosion  gas  bubble,  would  require  either  an 
unacceptably  coarse  mesh  or  an  unacceptably  large  number  of  finite  volumes.  This  program 
also  provides  a  constitutive  equation  suitable  for  modeling  of  explosive  charge  detonation, 
and  a  detonation  wave  front  algorithm  for  detonating  an  explosive  material. 
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This  processor  uses  the  control  volume  method  to  solve  the  basic  conservation 
equations  in  space.  The  integral  form  of  these  equations,  for  conservation  of  mass,  linear 
momentum,  and  total  energy,  are 

-w  J  J  J  volume  J  J  surface  *  ' 


f///        pudV=-f(        piKS'd§)*ff        TdS  (2-2) 
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-/// ,  ?edv=-H  *  f><*<®+n  „ aTd§         (2-3) 

QfJJJ  volume  •>  ->  surface  •>  J  surface  v  ' 

where  T  is  the  stress  tensor.  For  the  hydrodynamic  material  models  assumed  for  the 
numerical  analyses  of  this  investigation,  T  has  the  value  -p  along  the  main  diagonal,  and  zero 
everywhere  else.  Viscosity  and  thermal  conductivity  are  neglected  in  these  equations. 

These  equations  are  solved  for  each  finite  volume.  A  one  point  approximation  is  used 
to  solve  these  equations  (the  value  at  the  geometric  center  of  the  control  volume),  in 
conjunction  with  interpolated  velocities  and  pressures  at  the  faces  of  the  finite  volumes. 

A  first  order  "Donor-Acceptor"  scheme  is  used  for  material  transport.  Transported 
quantities  are  subtracted  from  donor  cells  and  added  to  the  acceptor  cells,  based  upon  the 
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donor  cell  values  and  the  velocity  at  the  common  face.  The  face  velocity  for  the  face 
connecting  finite  volumes  m  and  n  is 

and,  during  time  step  dt,  the  volume  transport  is 
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dV=afac;dSdt  (2-5) 

The  normal  component  of  the  face  velocity  determines  which  cell  is  a  donor  and  which  is  an 
acceptor  (i.e.  in  which  direction  flow  occur).  If,  for  example,  cell  m  is  the  donor,  the  mass, 
momentum,  and  energy  transport  from  element  m  to  n  is 

dM*Kdv=9jfa„'<$dt  (2-6) 

d(MU)=  Pjmdv=  9mam{aface-dS)dt  (2-7) 

*-P.V*r-P.«ai*ita,*fi*  (2-8) 

The  solution  in  time  is  computed  with  this  program  using  an  explicit  central  finite 
difference  method.  With  this  method,  a  new  time  step  size  is  calculated  after  each  forward 
time  step.  This  new  time  step  is  calculated  such  that  the  Courant  criterion  for  a  stable 
solution. 


dt=S—  (2-9) 

u+c 


is  satisfied,  where  S  is  a  factor  of  safety,  L  is  the  smallest  element  dimension,  c  is  the  acoustic 
wave  velocity  within  the  element,  and  u  is  the  partical  velocity  within  the  element.  This  is 
calculated  for  all  elements,  and  the  smallest  value  obtained  is  used  as  the  new  time  step. 

A  complete  time  step  algorithm  consists  of  several  phases.  The  momentum  change 
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due  to  the  impulse  (from  the  face  pressures)  is  calculated  and  used  to  update  the  velocity  in 
each  finite  volume.  Similarly,  the  work  done  by  the  pressure  on  each  finite  volume  face  is 
calculated  and  used  to  update  the  total  energy  for  each  finite  volume.  Mass,  momentum  and 
total  energy  transport  across  all  faces  are  computed,  and  used  to  update  these  values  in  all 
finite  volumes.  The  density  and  velocity  in  each  finite  volume  are  then  updated  (e.g.  density 
equals  new  mass  divided  by  fixed  volume).  The  specific  internal  energy  in  each  finite  volume 
is  computed  based  upon  the  new  specific  total  energy  and  the  new  specific  kinetic  energy. 
And  the  pressure  in  each  finite  volume  is  updated  using  the  new  density  and  specific  internal 
energy  in  the  state  equation. 

To  define  a  multimaterial  Eulerian  model  with  this  computer  program,  several  inputs 
are  required.  In  addition  to  defining  the  geometry  of  the  mesh,  constitutive  equation  data 
must  be  input  for  the  different  materials  to  be  used.  Initial  conditions  must  be  assigned  to  all 
of  the  finite  volumes  in  the  model;  for  a  multimaterial  Eulerian  problem  these  would  include 
specification  of  what  materials  are  initially  in  what  finite  volumes,  and  what  the  initial  density 
and  specific  internal  energy  is  in  the  different  regions  within  the  model.  Any  boundary 
conditions  which  are  not  to  be  left  as  the  default  "no  flow"  boundary  condition  must  be 
specified.  And,  if  the  problem  will  use  an  explosive  material,  one  or  more  detonation  points, 
detonation  velocities,  and  detonation  times  must  be  specified.  Finally,  the  ending  time  for  the 
analysis  must  be  specified,  and  the  program  must  be  told  precisely  which  output  quantities  are 
desired  at  what  frequency. 
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HI.  BASIC  EXPLOSION  GAS  BUBBLE  BEHAVIOR 

The  effect  of  fluid  compressibility  and  gas  internal  energy  on  the  behavior  of  a  free- 
field  explosion  gas  bubble  is  investigated  in  this  chapter.  As  discussed  in  the  introduction,  the 
basic  characteristics  of  an  explosion  gas  bubble  are  a  rapid  expansion  to  a  much  larger  volume 
than  the  initial  volume  of  the  charge,  a  prolonged  interval  of  time  in  this  expanded  state,  and 
a  rapid  collapse  back  to  a  small  volume.  This  oscillation  cycle  then  repeats,  but  with  a 
decreasing  period  between  collapses  and  a  decrease  in  the  maximum  volume  to  which  the 
bubble  expands. 

In  the  first   section  in  this  chapter,   the  basic  theoretical   equations  for  the 
noncompressive  motion  of  a  free-field  explosion  gas  bubble  are  developed.    The  second 
section  describes  the  details  of  the  numerical  analysis  models  used  in  this  investigation. 
Results  from  these  analyses  are  presented  in  the  last  section. 
A.         THEORETICAL  DEVELOPMENT 

The  classical  expression  for  the  motion  of  an  explosion  gas  bubble  in  the  absence  of 
a  nearby  boundary  surface  neglects  buoyancy  resulting  from  pressure  variation  around  the 
bubble  and  compressibility  of  the  surrounding  fluid.  With  these  restrictions,  the  flow  of  fluid 
outside  the  bubble  is  radial,  and  conservation  of  energy  gives 


/  .  \ 


yPo*3 


*K|   +^P(Ri+E(R)=Y  (3-1) 
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where  R  is  the  radius  of  the  bubble,  p0  is  the  density  in  the  surrounding  fluid,  PQ  is  the 
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hydrostatic  pressure  in  the  surrounding  fluid,  E(R)  is  the  internal  energy  of  the  gas  bubble 
when  it  has  radius  R,  and  Y  is  the  total  energy  (a  constant)  (Cole,  1948,  p. 273). 

In  this  expression,  the  first  term  represents  the  kinetic  energy  of  the  flow  in  the 
surrounding  fluid,  and  the  second  term  is  the  work  done  against  hydrostatic  pressure  in 
expanding  the  bubble  to  radius  R.  By  neglecting  the  internal  energy  of  the  gas,  which  has  a 
relatively  small  value  over  much  of  the  oscillation  cycle  of  an  explosion  gas  bubble,  the  total 
energy  can  be  expressed  as 

^y^l  (3-2) 

where  Rmax  is  the  maximum  radius  of  the  bubble  (Cole,  1948,  p. 274-275).  This  expression 
simply  reflects  the  fact  that  the  gas  bubble  must  have  a  maximum  radius  when  dR/dt  is  zero. 
Using  equation  (3-2)  in  equation  (3-1),  with  E(R)  taken  as  zero,  equation  (3-1)  can  be 
separated  and  integrated  to  give 


t=(3p0/2P0)1/2  h(RmJaf-  \\mda  (3.3) 


where  Rq  is  the  initial  radius  of  the  gas  bubble.  By  taking  R0  as  zero  (the  initial  radius  of  the 
bubble  is  much  smaller  than  the  maximum  radius)  and  transforming  and  integrating  equation 
(3-3),  Willis  (1941)  arrived  at  the  result 


T=l.Z3R 

max 


(3-4) 
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where  T  is  the  first  oscillation  period  of  the  bubble.  Substituting  equation  (3-2)  into  (3-4), 
and  using  the  fact  that  the  total  energy  Y  is  proportional  to  the  weight  W  of  the  explosive 
charge  gives 


T=Kr-^Te  (3-5) 


where  KT  is  a  constant  particular  to  a  given  type  of  explosive.   Equation  (3-2)  can  be  re- 
arranged to  give 

Rmax=KR—^  (3-6) 

ro 

where  KR  is  a  constant  particular  to  a  given  type  of  explosive.  Equations  (3-5)  and  (3-6)  are 
the  classical  "semi-empirical "  expressions  for  the  scaling  relationships  for  the  first  oscillation 
period  and  maximum  radius  of  an  explosion  gas  bubble. 
B.         NUMERICAL  ANALYSIS  MODELS 

Two  different  numerical  analyses  were  conducted  for  deep  explosion  gas  bubbles  in 
the  absence  of  a  nearby  boundary  surface.  The  first  analysis  was  designed  to  correspond  to 
one  of  a  series  of  experiments  conducted  by  Swift  and  Decius  (1950)  at  the  Woods  Hole 
Oceanographic  Institution  during  and  shortly  after  the  Second  World  War.  These  researchers 
conducted  a  large  number  of  experiments  to  determine  constants  for  the  semi-empirical 
equations  (3-5)  and  (3-6)  above  for  different  types  of  explosives.  The  tests  were  conducted 
at  a  deep  depth  in  deep  water,  to  minimize  the  effects  of  buoyancy  and  boundary  surfaces. 

The  problem  geometry  for  this  analysis  is  shown  in  Figure  3-1 .  A  very  small  (0.299 
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Figure  3-1.  Geometry  of  Problem  for  Free-Field  Bubble  Analysis 


kg  TNT  equivalent)  charge  was  detonated  at  a  depth  of  178.6  m  in  seawater.  The  0.299  kg 
TNT  equivalent  charge  weight  was  calculated  by  Swift  and  Decius  to  account  for  the  extra 
energy  of  the  detonator  and  booster  material  used  above  the  energy  of  the  TNT  in  the  main 
charge.  This  particular  shot  (shot  G72F)  was  chosen  from  all  of  the  shots  conducted  because 
the  bubble  oscillation  period  and  maximum  bubble  radius  from  this  experiment  were  closest 
to  the  mean  values  from  all  of  the  experiments.  Although  the  experiments  used  cylindrical 
charges  with  a  height  to  diameter  ratio  slightly  larger  than  one,  a  spherical  charge  was 
modeled  for  this  analysis. 

In  this  experiment  the  maximum  radius  of  the  bubble  was  much  less  than  the  depth  at 
which  it  was  located.  Thus,  in  the  numerical  analysis  model  for  this  problem,  gravity  was 
neglected  and  the  pressure  was  assumed  to  be  uniform  in  all  directions  from  the  charge.  With 
this  assumption,  the  problem  is  spherically  symmetric.  The  seawater  surrounding  the  charge 
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was  treated  as  an  invicid,  irrotational  fluid,  and  heat  and  mass  transfer  were  assumed  to  be 
negligible  over  the  time  frame  of  the  analysis.  Herring  (1950)  has  shown  that  the  transfer  of 
heat  over  one  bubble  oscillation  cycle  through  conduction  at  the  surface  of  the  bubble 
constitutes  a  negligible  fraction  of  the  total  energy  of  the  explosion. 

Since  this  problem  as  modeled  has  spherical  symmetry,  a  one-dimensional  model  is 
appropriate.  In  order  to  model  this  one-dimensional  problem  with  the  three-dimensional 
analysis  code  MSC/DYTRAN,  a  tall,  thin  pyramid  shaped  region  of  fluid  was  used.  This 
accounts  for  the  increasing  volume  as  one  moves  away  from  the  charge.  Boundary  conditions 
on  the  sides  of  this  fluid  volume  were  made  "no  flow"  conditions,  as  flow  across  the 
boundaries  is  precluded  by  the  flow  in  adjacent  fluid  volumes.  The  analysis  model  used  for 
this  problem  is  shown  in  Figure  3-2. 

At  the  apex  of  the  pyramid  shaped  volume,  a  rectangular  parallelpiped  was  used,  with 
the  equivalent  volume  a  pyramid  shaped  volume  would  have.  This  avoids  the  use  of  a  five 
sided  pyramid  shaped  volume  at  the  point,  which  is  not  a  standard  finite  element  shape. 
Rectangular  parallelpiped  shaped  finite  volumes  were  used  to  model  the  charge,  and  non- 
rectangular  hexahedron  finite  volumes  were  used  to  model  the  remaining  volume.  Because 
the  primary  interest  in  this  analysis  was  in  the  bubble  behavior  rather  than  the  primary  shock 
wave,  only  three  finite  volumes  were  used  to  model  the  charge. 

In  order  to  prevent  reflection  from  the  far  end  of  the  finite  volume  of  fluid  modeled 
from  affecting  the  results  during  the  duration  of  the  analysis,  a  very  large  volume  of 
surrounding  fluid  was  modeled.  A  total  of  999  finite  volumes  were  used  for  this  analysis 
model,  with  the  radial  length  of  individual  elements  increasing  with  increasing  distance  from 
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COMPLETE  MODEL 


REGION  NEAR  CHARGE 


Figure  3-2.  Numerical  Analysis  Model  for  0.299  kg  TNT  Explosive  Charge 
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the  charge. 

The  TNT  for  this  problem  was  modeled  using  a  JWL  state  equation,  with  state 
equation  parameters  taken  from  the  Lawrence  Livermore  National  Laboratory  Explosives 
Handbook  (Dobratz,  1981).  With  this  state  equation,  the  pressure  in  the  "burned  fraction" 
of  a  charge  is  related  to  the  specific  internal  energy  and  density  by 


p=A 
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(3-7) 


where  r|=p/po,  p0  is  the  initial  density,  E  is  the  specific  internal  energy  (per  unit  mass),  and  A, 

B,  co,  Rlt  and  R2  are  constant  parameters  for  the  explosive.  The  "burned  fraction"  is  just  that 

portion  of  the  explosive  contained  within  a  spherical  detonation  front  traveling  outward  with 

detonation  velocity  d,  which  is  another  parameter  characterizing  the  explosive.  For  TNT 

compressed  to  an  initial  density  of  p0  of  1 .630  g/cm3,  the  JWL  state  equation  parameters  are 

(Dobratz,  1981): 

£  =  4.29xl06J/kg 
A  =  3.712  xlO11  Pa 
5  =  3.231  xlO9  Pa 
co  =  0.30 
#,  =  4.15 
R2  =  0.95 
</=6930m/s 

In  order  to  model  the  seawater  in  which  this  experiment  was  conducted,  a  polynomial 
state  equation  was  used.  This  state  equation  relates  the  pressure  in  the  fluid  to  its  density  p 
and  specific  internal  energy  E  per  unit  mass  by 
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p=alii+a2ii2+a3n:i+(b0+blii+b2\i2)i>GE  (3-8) 


where 


P-Po 
H= (3-9) 


is  the  condensation  and  alt  a2,  a3,  b0,  bx,  and  b2  are  constants  for  the  fluid.  Constant  for  this 
state  equation  for  seawater,  for  use  at  condensation  values  u  of  up  to  0.8,  were  determined 
by  fitting  available  Hugonoit  data  from  the  literature  to  this  state  equation  form,  and  adjusting 
the  bulk  modulus  and  density  to  the  accepted  values  for  seawater.  The  resultant  state 
equation  parameters  were  (Chisum  and  Shin,  1995) 

^9 


a,  =  2.306x10*  Pa 
h 


a,  =  8.432  xlO9  Pa 


a3  =  8.014  xlO9  Pa 
b0  =  0.4934 
6,  =  1.3937 
b2  =  0.0000 
p0=1025kg/m3 

The  initial  conditions  given  to  the  seawater  in  this  problem  were  an  initial 

condensation  value  of  zero  (initial  density  of  1025  kg/m3),  and  an  initial  specific  internal 

energy  of  3750.4  J/kg.  This  initial  specific  internal  energy  was  determined  from  equation  (3- 

8);  it  represents  the  specific  internal  energy  necessary  to  give  the  seawater  an  initial  pressure 

equal  to  the  hydrostatic  plus  atmospheric  pressure  for  this  problem.  It  was  found  necessary 

to  use  the  initial  specific  internal  energy  rather  than  the  condensation  to  set  the  initial 

pressure,  because  seawater  is  so  incompressible  that  there  is  only  a  minuscule  density  change 

at  this  depth,  and  significant  round  off  errors  would  be  introduced  by  setting  the  initial 
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pressure  with  the  condensation. 

The  second  model  geometry  analyzed  was  similar  to  Figure  3-1,  except  that  an  initial 
charge  depth  of  1000  m  and  an  initial  TNT  charge  of  3.5  kg  was  used.  This  corresponds  to 
an  initial  charge  radius  of  8.0  cm.  The  objective  of  this  analysis  was  to  examine  the  flow 
characteristics  in  the  surrounding  fluid  over  one  oscillation  of  the  bubble.  For  this  reason, 
uniform  element  spacing  in  the  radial  direction  was  used.  A  pyramid  shaped  volume  of  fluid 
was  again  used.  In  this  case,  the  fluid  was  only  modeled  out  to  a  distance  of  16  m,  and  a 
simplified  "non-reflecting"  boundary  condition  was  placed  on  the  last  element.  This  boundary 
condition  was  that  the  flow  across  the  face  of  the  element  at  the  boundary  would  have  the 
same  values  as  in  that  element. 

The  radial  element  spacing  used  was  0.125  cm,  and  a  total  of  12800  finite  volumes 
were  used.  It  was  not  necessary  to  use  this  many  elements  to  get  a  "grid  independent" 
solution,  but  the  use  of  a  fine  mesh  better  illustrates  the  characteristics  of  the  shock  wave. 
The  computer  program  used  utilizes  artificial  bulk  viscosity  to  control  oscillations  behind 
shock  waves,  and  this  smears  the  shock  front  over  a  number  of  elements.  Mass,  energy,  and 
momentum  are  still  conserved,  but  the  shock  wave  doesn't  look  as  steep.  Thus,  using  more 
elements  makes  the  shock  wave  "look"  more  like  a  shock  wave,  since  it  is  smeared  over  a 
smaller  distance. 
C.         NUMERICAL  ANALYSIS  RESULTS 

As  stated  in  the  last  section,  the  first  of  the  problems  analyzed  corresponded  to  an 
experiment  conducted  by  Swift  and  Decius.  Figure  3-3  shows  the  predicted  radius  time 
history  from  this  analysis,  and  the  first  maximum  bubble  radius  and  period  measured  from  the 
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Figure  3-3.  Radius  Time  History  for  0.299  kg  TNT  Charge  Detonated  at  178.6  m 
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experiment.  This  plot  shows  excellent  agreement  between  the  analysis  results  and  the 
experimental  measurement.  Also  shown  on  this  plot,  for  comparison  purposes,  is  the  radius 
time  history  curve  generated  from  step-by-step  numerical  integration  of  equation  (3-3),  which 
neglects  both  the  internal  energy  of  the  gas  bubble  and  fluid  compressibility.  This  curve  is 
seen  to  agree  quite  well  with  the  numerical  analysis  curve  up  until  the  bubble  reaches 
maximum  volume,  but  to  begin  diverging  from  there. 

It  is  seen  in  Figure  3-3  that  the  numerical  analysis  predicts  the  bubble  radius  and 
oscillation  period  to  decrease  in  subsequent  oscillations.  This  is  due  to  the  radiation  of  a 
bubble  pulse  when  the  bubble  collapses.  This  effect  cannot  be  generated  without  including 
fluid  compressibility. 

Figure  3-4  shows  the  pressure,  impulse,  and  fluid  particle  velocity  time  histories,  at 
a  point  located  two  maximum  bubble  radii  from  the  center  of  the  charge.  This  figure 
illustrates  that,  although  the  peak  pressure  in  the  first  bubble  pulse  is  much  less  than  the  peak 
shock  wave  pressure,  the  impulse  generated  by  the  bubble  pulse  is  still  large  because  of  its 
longer  duration. 

In  their  experiment,  Swift  and  Decius  also  measured  the  second  period  and  second 
maximum  radius  for  this  shot.  Their  measurements  indicated  that  the  second  maximum  radius 
was  1 1 .6  inches  (about  29.5  cm),  and  that  the  measured  time  for  the  second  bubble  minimum 
was  30.85  msec.  These  results  indicate  that  an  additional  energy  loss,  above  that  due  to 
acoustic  radiation,  occurred  in  this  experiment,  since  acoustic  radiation  was  accounted  for  in 
the  numerical  analysis  for  this  problem. 

One  possible  explanation  for  this  in  that  the  bubble  surface  was  unstable  near  the 
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Figure  3-4.  Numerically  Predicted  Pressure,  Impulse,  and  Velocity  Time  History  for 
0.299  kg  TNT  Charge  Detonated  at  178.6  m,  at  a  Point  Two  Maximum 
Bubble  Radii  From  the  Charge  Center 
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minimum  volume.  Hicks  (1970a)  has  noted  that  photographs  of  non-migrating  bubbles 
typically  show  numerous  needle-like  water  jets  protruding  into  the  bubble  surface  when  the 
bubble  is  near  minimum  radius,  and  that  this  "spray"  of  water  could  significantly  cool  the  hot 
bubble  gases. 

In  the  second  free-field  numerical  analysis,  the  flow  external  to  the  bubble  was 
examined.  Figure  3-5  shows  the  bubble  radius  time  history  for  this  analysis,  and  Figure  3-6 
shows  the  variation  of  the  internal  energy  of  the  bubble  as  a  function  of  time.  The  internal 
energy  is  seen  to  rapidly  decrease  from  its  initial  value  as  the  bubble  expands,  and  then  to 
increase  moderately  as  the  bubble  is  compressed  by  the  surrounding  fluid. 

Figures  3-7,  3-8,  and  3-9  show  the  pressure  in  the  fluid  at  1  msec  intervals.  Because 
the  shock  wave  is  radiating  spherically,  the  peak  pressure  of  the  shock  wave  decreases  rapidly 
with  increasing  distance.  The  pressure  behind  the  shock  front  drops  below  hydrostatic 
pressure,  and  this  low  pressure  region  extends  for  a  considerable  distance  from  the  bubble  by 
5  msec,  when  the  bubble  is  near  maximum  radius.  Also,  by  5  msec  there  is  a  second  sharp 
drop  in  the  pressure  behind  the  shock  wave,  and  after  5  msec  this  sharp  drop  extends  further 
and  further  into  the  fluid. 

This  is  caused  by  the  reversal  of  flow  in  a  portion  of  the  fluid.  The  fluid  immediately 
behind  the  shock  front  is  still  moving  outward,  but  the  fluid  near  the  bubble  surface  is  moving 
inward  as  the  bubble  contracts.  As  time  goes  on,  fluid  further  and  further  away  from  the 
bubble  begins  moving  inward.  By  the  time  the  bubble  reaches  minimum  volume,  the  pressure 
near  the  bubble  is  significantly  higher  than  even  the  primary  shock  wave.  As  the  bubble 
reaches  minimum  volume  and  begins  expanding  again,  a  moderately  low  pressure  but  long 
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Figure  3-5.  Predicted  Bubble  Radius  Time  History  for  3.5  kg  TNT  Charge  Detonated  at 
1000  m 
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Figure  3-6.  Predicted  Bubble  Internal  Energy  Time  History  for  3.5  kg  TNT  Charge 
Detonated  at  1000  m 
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Figure  3-7.  Pressure  Distribution  From  1  msec  Through  5  msec  After  Detonation  of  a 
3.5  kg  TNT  Charge  at  a  Depth  of  1000  m 
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Figure  3-8.  Pressure  Distribution  From  5  msec  Through  9  msec  After  Detonation  of  a 
3.5  kg  TNT  Charge  at  a  Depth  of  1000  m 
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Figure  3-9.  Pressure  Distribution  From  9  msec  Through  13  msec  After  Detonation  of  a 
3.5  kg  TNT  Charge  at  a  Depth  of  1000  m 
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duration  bubble  pulse  begins  traveling  outward  away  from  the  bubble. 

Table  3-1  compares  the  results  between  the  experiment  conducted  by  Swift  and 
Decius  and  the  current  numerical  analysis  predictions.  Extremely  good  agreement  is  seen  for 
the  first  maximum  bubble  radius  and  period,  while  the  second  maximum  bubble  radius  and 
period  are  appreciably  less.  The  solution  computed  using  the  noncompressive  theory 
neglecting  the  internal  energy  of  the  gas  bubble  was  8%  lower  than  the  actual  period.  Internal 
gas  energy  and  fluid  compressibility  thus  have  an  observable  effect  on  the  results. 


EXPERIMENTAL 
MEASURMENT 

NUMERICAL 
ANALYSIS 

% 
ERROR 

R(l) 

39.1  cm 

38.8  cm 

0.8 

T(l) 

17.85  msec 

17.77  msec 

0.4 

R(2) 

29.5  cm 

32.2  cm 

9 

T(2) 

13.00  msec 

15.85  msec 

22 

Table  3-1 .  Comparison  Between  Experimental  and  Numerical  Analysis  Results  For  First  and 
Second  Maximum  Bubble  Radius  and  First  and  Second  Bubble  Period. 
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IV.  EXPLOSION  GAS  BUBBLE  BEHAVIOR  NEAR  SIMPLE  BOUNDARIES 

A.         INTRODUCTION 

While  free-field  explosion  gas  bubbles  dynamic  behavior  has  been  discussed  in  the 
previous  chapter,  practical  problems  typically  involve  the  interaction  of  these  bubbles  with 
boundary  surfaces.  As  discussed  in  the  introduction,  additional  phenomena  are  observed 
when  an  explosion  gas  bubble  is  generated  near  a  boundary  surface.  These  include  a  change 
in  the  oscillation  period  of  the  bubble,  and  migration  of  the  bubble  from  its  initial  position  in 
the  fluid  medium. 

These  phenomena  will  affect  the  damage  caused  by  whipping  of  the  hull  girder  of  a 
naval  vessel.  The  time  interval  between  arrival  of  the  primary  shock  wave  and  the  first  bubble 
pulse  can  be  very  near  or  equal  to  the  fundamental  period  of  the  whipping  mode  of  vibration 
of  the  hull  girder,  significantly  increasing  the  whipping  mode  response  of  the  vessel.  Since 
the  bubble  pulse  is  generated  at  the  conclusion  of  the  contraction  phase  of  the  bubble 
oscillation,  the  period  of  this  oscillation  is  obviously  important.  The  migration  of  the  bubble 
towards  or  away  from  a  boundary  during  this  oscillation  is  also  important,  as  this  determines 
the  origination  point  of  the  bubble  pulse.  This  affects  both  the  arrival  time  of  the  bubble  pulse 
and  its  magnitude,  as  the  magnitude  decreases  with  increasing  distance. 

Since  the  available  analytical  models  are  known  to  not  give  quantitatively  accurate 
results  for  the  effects  of  boundary  surfaces  on  underwater  explosion  gas  bubbles,  a  number 
of  numerical  analyses  using  Eulerian  finite  volume  meshes  were  conducted.  These  included 
a  free-field  analysis  in  which  no  nearby  boundary  was  present,  a  series  of  analyses  in  which 
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Charge   (10.24  kg   TNT) 
Free-Field   Maximum  Radius 

P    (atm)    =    101325   Pa 
P   =    1025   kg/m3 
g    =    9.80665   m/s2     I 
h   =    1000  m 

P^  =   P   (atm)    +   pgh 


Boundary   (various   standoffs   shown) 


Figure  4-1.  Geometry  of  Problem  for  Bubble  Near  Boundary  Analyses 


a  plane  infinite  rigid  boundary  was  located  at  various  standoff  distances  from  the  charge,  and 
a  series  of  analyses  in  which  a  plane  infinite  free  (constant  pressure)  surface  was  located  at 
various  standoff  distances  from  the  charge. 


B. 


NUMERICAL  ANALYSIS  MODELS 


The  basic  problem  geometry  for  these  analyses  consists  of  a  10.24  kg  cylinder  of 
TNT,  with  a  height  to  diameter  ratio  of  one,  located  near  an  infinite  plane  boundary  at  a  depth 
of  1000  m,  as  shown  in  Figure  4-1.  Underwater  explosion  experiments  also  typically  use 
cylindrical  charges,  for  practical  reasons.  As  in  the  free-field  analyses  described  in  the 
previous  chapter,  the  choice  of  a  deep  charge  depth  serves  to  significantly  reduce  the  period 
of  oscillation  of  the  bubble,  and  hence  the  analysis  time. 

A  deep  charge  depth  also  serves  to  simplify  the  analyses,  as  at  a  deep  depth  the 
variation  in  hydrostatic  pressure  in  the  surrounding  fluid  for  a  given  depth  change  is 
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proportionately  much  smaller  than  at  a  shallow  depth.  The  hydrostatic  pressure  in  the  fluid 
near  a  deep  charge  can  thus  be  approximated  as  having  a  uniform  value,  the  hydrostatic 
pressure  at  the  depth  of  the  charge,  provided  that  the  overall  dimensions  of  the  explosion  gas 
bubble  is  much  smaller  than  the  depth  of  the  charge.  This  was  the  assumption  made  for  the 
analyses  described  in  this  chapter,  for  which  the  maximum  bubble  radius  was  less  than  0. 1% 
of  the  depth  of  the  charge. 

The  migration  due  to  gravity  can  be  separated  from  the  migration  caused  by  the 
presence  of  plane  boundary  surfaces,  for  cases  in  which  the  boundary  surfaces  have  no 
component  normal  to  the  direction  of  the  gravitational  acceleration.  Furthermore,  deep 
explosion  gas  bubbles  are  known  to  experience  little  vertical  migration  due  to  gravity  (Hicks, 
1970a).  Thus,  gravity  can  be  neglected  except  for  its  effect  on  the  hydrostatic  pressure  in  the 
surrounding  fluid.  This  permits  separation  of  the  effects  of  gravity  from  the  effects  of  a 
boundary  surface,  facilitating  investigation  of  the  influence  of  the  boundary  surface  on  the 
dynamic  behavior  of  the  bubble.  For  the  same  purpose,  Blake  and  Gibson  (1981)  have 
conducted  experiments  using  spark  generated  vapor  bubbles  with  the  entire  experimental 
apparatus  in  free  fall. 

For  these  analyses,  the  cylindrical  charge  is  assumed  to  have  its  axis  normal  to  the 
plane  boundary  surfaces.  With  this  orientation,  the  problems  analyzed  are  axially  symmetric, 
with  the  symmetry  axis  being  the  axis  of  the  charge.  To  analyze  these  axially  symmetric 
geometries  using  the  three-dimensional  code  MSC/DYTRAN,  wedge  shaped  volumes 
comprising  two  degrees  of  arc  were  used. 

Figure  4-2  illustrates  the  overall  shape  of  the  volume  of  fluid  modeled,  for  the  free- 
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Figure  4-2.  Overall  Geometry  of  Analysis  Model  for  Free-Field  Analysis 
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field  analysis  in  which  no  boundary  surface  is  present.  The  vertical  line  oriented  in  the  y 
direction  in  this  model  is  the  axisymmetric  symmetry  axis  of  the  problem.  The  boundary 
conditions  imposed  on  the  two  semi-circular  areas,  lying  in  the  planes  oriented  at  +/-  1  degree 
about  the  axisymmetric  symmetry  axis  from  the  x-y  plane,  are  "no-flow"  (rigid  wall)  boundary 
conditions,  as  flow  across  these  boundaries  is  precluded  by  the  flow  in  adjacent  wedge 
segments  (not  modeled  or  shown). 

As  discussed  in  the  section  in  the  previous  chapter  describing  free-field  analysis 
results,  the  dynamic  behavior  of  underwater  explosion  gas  bubbles  is  influenced  by  the 
conditions  in  the  surrounding  fluid  adjacent  to  and  at  some  distance  from  the  bubble.  Thus, 
a  large  volume  of  surrounding  fluid  was  modeled  in  these  analyses.  This  acts  both  to 
incorporate  the  necessary  fluid  and  to  avoid  having  to  specify  particular  boundary  conditions 
at  the  remote  boundary  of  the  finite  volume  of  modeled  fluid.  This  remote  boundary  was 
simply  made  a  "no-flow"  (rigid  wall)  boundary,  with  the  distance  to  the  boundary  being  made 
so  great  that  shock  wave  reflection  from  it  could  not  affect  the  analyses  during  the  time  frame 
of  the  analyses.  For  the  free-field  model  volume  shown  in  Figure  4-2,  this  remote  boundary 
is  the  curved  area  connecting  the  semi-circles  rotated  +/-  1  degree  about  thej-axis  from  the 
x-y  plane.  These  semi-circles  are  located  at  a  distance  of  400  m  from  the  detonation  point; 
as  discussed  below,  an  increased  element  size  is  used  far  away  from  the  center  of  the  charge, 
so  that  relatively  few  extra  elements  are  needed  to  model  fluid  very  remote  from  the  center 
of  the  charge. 

The  finite  volume  meshes  used  for  the  analyses  described  in  this  chapter  were 
developed  to  meet  several  goals.  It  was  important  that  the  results  obtained  be  independent 
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of  the  mesh  used  to  obtain  them,  i.e.  that  the  solution  was  "grid  independent."  This  goal 
could  be  met  simply  by  using  a  large  number  of  equal  sized  elements  (finite  volumes)  for  the 
analyses.  However,  if  carried  to  extremes  this  practice  is  in  direct  conflict  with  other  goals 
that  a  practical  model  must  meet;  that  is,  that  the  model  must  be  able  to  run  using  a 
reasonable  amount  of  computer  resources  (such  as  memory  and  storage  space  for  results), 
and  that  the  model  must  run  to  completion  in  a  reasonable  amount  of  time.  It  was  also 
considered  important  that  the  results  be  independent  of  the  shape  of  the  elements  used  in  the 
area  near  the  charge,  even  if  the  bubble  migrated  during  the  analyses.  A  considerable  amount 
of  effort  was  expended  in  order  to  meet  all  of  these  goals  for  the  analyses  described  in  this 
chapter. 

During  the  development  of  meshes  for  these  analysis,  it  was  eventually  determined 
that  a  mesh  consisting  of  three  separate  regions  could  meet  these  goals.  From  a  two- 
dimensional  perspective,  the  central  region  of  these  meshes  consists  of  a  rectangular  region 
with  square  elements.  This  configuration  allows  elimination  of  element  shape  as  a  factor, 
since  if  this  element  geometry  can  still  produce  spherically  shaped  bubbles  and  the  bubbles 
remain  within  this  region  during  the  analyses,  then  the  shape  of  the  element  was  not  important 
in  determining  the  results.  The  second  region  is  a  transition  from  this  square  element  region 
to  a  spherically  diverging  type  of  region,  which  is  necessary  if  the  total  number  of  elements 
is  to  be  kept  to  a  reasonable  value.  The  third  region  is  the  spherically  diverging  type  region, 
in  which  the  dimensions  of  the  element  increase  with  increasing  distance  from  the  center  of 
the  charge,  which,  again,  serves  to  limit  the  number  of  elements  to  a  reasonable  value. 

The  use  of  this  spherically  diverging  region  is  supported  and  suggested  by  the  results 
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from  the  previous  chapter,  where  it  is  seen  that  characteristics  of  the  flow  in  the  surrounding 
fluid  are  changing  much  slower  in  the  remote  fluid  than  in  the  fluid  near  the  charge.  As  the 
purpose  of  the  analyses  described  in  this  chapter  is  to  determine  the  dynamic  behavior  of  the 
bubble  rather  than  the  primary  shock  wave,  it  does  no  harm  if  the  primary  shock  wave  is 
spread  over  a  larger  region.  And  since  a  reasonable  number  of  elements  is  still  used  in  the 
circumferential  direction  (in  the  free-field  model,  34  elements  over  45  degrees,  i.e.  about  1.3 
degrees  per  element),  this  element  geometry  should  still  be  valid  for  the  analyses  in  which  a 
boundary  surface  is  present,  and  there  is  some  circumferential  change  in  the  flow  field. 

The  dimensions  of  the  rectangular  (square  meshed)  region  used  was  several  times  the 
maximum  bubble  radius  of  the  charge.  This  was  found  to  be  necessary  to  eliminate  mesh 
reflection  affects  on  the  bubble  radius-time  history  for  the  free-field  analysis.  This  meshing 
in  sufficient  to  capture  the  important  effects  which  occur  within  a  short  distance  of  the  charge, 
as  seen  in  the  previous  chapter. 

The  mesh  transition  region  provides  a  transition  from  the  square  meshed  central  region 
to  the  spherically  diverging  outer  region,  and  reduces  the  number  of  elements  in  the 
circumferential  direction  per  unit  circumferential  distance.  This  reduction  is  important  in 
limiting  the  total  number  of  elements  needed  for  the  analyses.  Figure  4-3  shows  the  geometry 
of  the  central  region  and  transition  region  for  the  free-field  analysis  model;  this  is  just  an 
expanded  view  of  Figure  4-2,  in  the  region  in  which  the  charge  is  located. 

To  insure  a  grid  independent  solution,  free-field  analyses  were  conducted  for  which 
the  square  element  dimensions  in  the  central  region  were  2.5  cm  and  3.333  cm  (with  slight 
differences  in  the  gross  dimensions  of  the  central  region  (247.5  cm  and  250  cm  width)),  and 
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Figure  4-3.  Geometry  of  Analysis  Model  for  Free-Field  Analysis  in  Area  Near  Charge 
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the  bubble  radius  versus  time  curves  were  found  to  lie  on  top  of  each  other.  It  was  concluded 
from  these  results  that  the  2.5  cm  square  dimension  solution  was  a  grid  independent  solution 
for  the  free-field  case.  It  was  then  assumed  that  this  meshing  would  be  adequate  for  analyses 
involving  a  nearby  boundary  surface,  since  the  meshing  for  these  analyses  was  kept  similar 
to  the  free-field  meshing. 

Figure  4-4  shows  the  gross  finite  volume  analysis  mesh  from  a  two-dimensional 
perspective,  for  the  free-field  analysis.  Figure  4-5  shows  the  central  rectangular  region  of  this 
mesh,  and  the  transition  from  this  central  region  to  the  outer  (spherically  diverging)  region. 
The  total  number  of  finite  volumes  (elements)  used  for  this  analysis  is  47340.  The  central 
(square  meshed)  region  uses  19602  elements,  the  transition  region  uses  338  elements,  and  the 
outer  region  uses  the  remaining  27200  elements.  In  this  outer  region,  the  mesh  is  designed 
such  that  the  dimensions  of  individual  elements  remain  nearly  equal  as  the  mesh  diverges  from 
the  central  region. 

The  maximum  bubble  radius  predicted  by  the  free-field  analysis  was  about  70.65  cm, 
and  the  first  bubble  period  from  this  analysis  was  about  14.47  msec.  These  values  were  used 
to  assist  in  quantification  of  the  results  for  analyses  which  included  a  boundary  surface.  For 
analyses  involving  a  boundary  surface,  lengths  are  given  in  terms  of  the  non-dimensional 
standoff  distance  h*,  where  h*  is  the  standoff  distance  from  the  center  of  the  charge  to  the 
nearest  point  on  the  boundary  in  units  of  maximum  free-field  gas  bubble  radii.  Times  are 
given  in  terms  of  the  non-dimensional  period  T*,  where  T*  is  the  time  in  units  of  bubble  free- 
field  first  oscillation  periods. 

The  free-field  axisymmetric  model  was  used  as  the  starting  point  for  other  analysis 
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Figure  4-4.  Overall  Finite  Volume  Analysis  Mesh  for  Free-Field  Analysis 
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Figure  4-5.  Finite  Volume  Analysis  Mesh  for  Free-Field  Analysis  in  Area  Near  Charge 
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models  involving  a  rigid  wall  boundary  and  a  constant  pressure  surface.  Analyses  were 
conducted  at  standoff  distances  of  h*  ranging  from  1.062  to  4.034,  for  both  types  of 
boundaries.  These  h*  values  correspond  to  standoff  distances  from  the  center  of  the  charge 
to  the  boundary  ranging  from  75  cm  to  285  cm,  respectively. 

Table  4-1  summarizes  the  characteristics  of  the  models  used  in  these  analyses.  The 
total  number  of  finite  volumes  (cells)  used  in  these  models  ranged  from  28680  (for  the  h* 
equal  to  1.062  analyses)  to  42708  (for  the  h*  equal  to  4.034  analyses).  The  models  for  these 
analyses  were  similar  to  the  free-field  model  described  above,  but  with  additional  elements 
present  below  what  was  an  additional  symmetry  plane  for  the  free-field  model  (bisecting  the 
axis  of  the  charge). 


Standoff  Distance  (cm) 

ft* 

Number  of  Finite  Volumes 

75.0 

1.062 

28680 

82.5 

1.168 

29181 

97.5 

1.380 

30183 

120.0 

1.698 

31686 

142.5 

2.017 

33189 

187.5 

2.654 

36195 

285.0 

4.034 

42708 

(Free-Field) 

(Infinity) 

47340 

Table  4- 1 .   Summary  of  Characteristics  of  Finite  Volume  Analysis  Models 


Figure  4-6  shows  the  geometry  of  a  typical  model  used  in  these  analyses  from  a  three- 
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Figure  4-6.  Overall  Model  Geometry  for  h*=1.380  Analyses 
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dimensional  perspective.  This  figure  is  for  the  analyses  in  which  h*  has  a  value  of  1.380. 
Figure  4-7  shows  an  expanded  view  of  the  geometry  of  this  model,  near  the  central  region. 
In  these  figures,  the  lower  of  the  triangular  shaped  areas  parallel  to  the  z-x  plane  is  the  rigid 
wall  or  constant  pressure  boundary.  The  upper  (interior)  triangular  shaped  area  bisects  the 
center  of  the  charge,  and  was  the  symmetry  plane  for  the  free-field  analysis. 

The  overall  finite  volume  mesh  used  for  the  model  geometry  shown  in  Figures  4-6  and 
4-7  is  shown  in  Figure  4-8,  from  a  two-dimensional  perspective.  Figure  4-9  shows  an 
expanded  two-dimensional  view  of  this  model  in  the  region  near  the  charge.  As  in  the  free- 
field  model,  this  model  has  a  square  meshed  central  region,  which  for  this  model  is  a 
rectangular  wedge  247.5  cm  wide  by  345  cm  high.  This  is  bordered  on  two  sides  by  a 
transition  region  connecting  this  square  meshed  area  with  the  more  coarsely  meshed  outer 
region. 

The  lower  boundary  parallel  to  the  z-x  plane  in  Figures  4-6  through  4-9  is  made  either 
a  rigid  boundary  or  a  constant  pressure  boundary,  by  imposition  of  boundary  conditions  on 
the  faces  of  the  elements  which  lie  on  this  boundary.  For  the  rigid  wall  analysis,  the 
appropriate  boundary  condition  is  a  "no-flow"  condition.  For  the  constant  pressure  boundary 
analysis,  the  pressure  on  the  element  faces  which  constitute  this  boundary  is  made  equal  to 
the  initial  hydrostatic  pressure  in  the  surrounding  fluid,  for  all  time.  The  other  flow 
parameters  (velocity,  density,  and  specific  internal  energy)  at  the  boundary  then  take  on  the 
values  in  the  finite  volumes  in  to  or  out  of  which  material  is  flowing.  This  constant  pressure 
boundary  is  only  an  approximation  of  an  actual  air- water  interface,  as  it  does  not  account  for 
changes  in  the  shape  of  the  boundary  (which  is  assumed  to  remain  planar). 
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Figure  4-7.  Model  Geometry  for  h*=l  .380  Analyses  in  Area  Near  Charge 
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Figure  4-8.  Overall  Finite  Volume  Analysis  Mesh  for  h*=l  .380  Analyses 
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The  finite  volume  models  for  the  other  analyses  conducted  at  various  other  standoff 
distances  are  very  similar  to  the  model  shown  in  Figures  4-6  through  4-9.  These  models  vary 
only  in  the  amount  of  fluid  added  between  the  plane  bisecting  the  axis  of  the  charge  and  the 
boundary  plane. 

The  10.24  kg  cylindrical  TNT  charge  used  in  the  simulations  discussed  in  this  chapter 
is  modeled  as  a  20  cm  high,  20  cm  diameter  charge  with  an  initial  density  of  1.630  g/cm3. 
Properties  for  this  charge  were  modeled  using  a  JWL  form  state  equation,  with  state  equation 
parameters  taken  from  the  Lawrence  Livermore  National  Laboratory  Explosives  Handbook 
(Dobratz,  1981).  This  is  the  same  state  equation  described  in  Chapter  III,  and  the  same 
parameters  given  there  for  this  state  equation  were  used  for  the  analyses  described  in  this 
chapter.  In  Figures  4-8  and  4-9,  all  of  the  elements  initially  contain  seawater,  with  the 
exception  of  the  32  elements  initially  containing  the  TNT  charge.  This  charge  is  shown  as 
shaded  elements  in  Figure  4-10.  This  TNT  is  detonated  using  a  spherical  detonation  wave 
traveling  outward  from  the  center  of  the  charge  (midway  up  the  charge  on  the  axisymmetric 
symmetry  axis)  at  a  constant  velocity  of  6930  m/s. 

The  seawater  in  all  of  the  remaining  elements  was  modeled  using  the  same  polynomial 
state  equation  and  state  equation  parameters  described  in  Chapter  m.  Initial  conditions  in  this 
seawater  were  an  initial  density  of  1.025  g/cm3,  and  an  initial  pressure  of  10.153  MPa 
(0. 10153  Kbar),  which  is  the  total  hydrostatic  pressure  (atmospheric  plus  seawater  head).  As 
discussed  in  Chapter  III,  this  initial  pressure  was  set  by  specifying  a  nonzero  initial  specific 
internal  energy  (20076  J/Kg)  to  avoid  round  off  errors  (even  at  a  depth  of  1000  m,  there  is 
only  a  minuscule  density  increase).  Using  a  constant  volume  specific  heat  value  of  about  1 
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Figure  4-10.  Locate  of  Charge  in  Finite  Volume  Analysis  Mesh  for  h*=1.380  Analyses 
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cal/g  C,  this  works  out  to  a  temperature  change  of  only  7.8  C,  for  which  the  results  should 
not  be  affected. 

Several  simplifying  assumptions  were  made  in  these  numerical  analyses.  As  in  the 
simulation  described  in  the  last  chapter,  the  seawater  was  assumed  to  behave  as  an  invicid  and 
irrotational  fluid.  Heat  and  mass  transfer  between  the  seawater  and  the  explosion  gas  bubble 
were  also  assumed  to  be  negligible  over  the  time  frame  of  the  analyses.  Herring  (1950)  has 
shown  the  peak  stresses  due  to  viscosity  are  negligible  in  comparison  with  hydrostatic 
pressure,  even  using  Taylor's  (1943)  velocity  field  (which  is  known  to  significantly 
overestimate  the  bubble  migration  velocity),  and  that  the  transfer  of  heat  over  one  bubble 
oscillation  cycle  through  conduction  at  the  surface  of  the  bubble  constitutes  a  negligible 
fraction  of  the  total  energy  of  the  explosion. 
C.         SIMULATION  RESULTS 

1.  The  Effect  of  Boundaries  on  Bubble  Period  in  a  Compressible  Fluid 

As  discussed  in  the  introduction,  the  period  of  the  first  bubble  oscillation  and  the 
displacement  of  the  bubble  during  this  time  interval  are  important  factors  in  determining  the 
amount  of  damage  that  the  bubble  can  cause.  Figure  4-1 1  shows  the  volume  equivalent 
spherical  radius  (the  radius  of  a  spherical  bubble  with  the  same  volume)  time  history,  for 
analyses  in  which  a  rigid  and  a  constant  pressure  boundary  were  located  at  a  standoff  distance 
of  1.380  maximum  free-field  radii.  The  use  of  the  volume  equivalent  spherical  radius  as  the 
ordinate  scale  in  this  figure  is  not  meant  to  imply  that  bubbles  were  actually  spherical  when 
a  boundary  was  present.  This  ordinate  scale  is  simply  (3  V/47i)1/3,  and  was  chosen  merely  for 
convenience  in  comparing  results  from  the  various  analyses  with  the  free-field  analysis  results. 
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Figure  4-11.  Bubble  Radius  Time  History  for  Free-Field  and  h*=1.380  Analyses 
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The  free-field  radius-time  history  is  also  shown  in  this  graph  for  comparison  purposes.  The 
period  of  the  bubble  is  obviously  influenced  appreciably  by  the  presence  of  the  nearby 
boundary,  as  is  the  maximum  volume  of  the  bubble.  Similar  bubble  volume  equivalent  radius 
time  history  data  was  generated  for  the  analyses  conducted  at  the  other  initial  standoff 
distances,  and  used  to  determine  the  individual  bubble  oscillation  periods. 

By  compiling  the  bubble  period  results  from  the  numerical  analyses  conducted  at 
different  standoff  distances,  the  effect  of  the  standoff  distance  to  each  type  of  boundary  can 
be  assessed.  This  is  done  in  Figure  4-12,  in  which  the  non-dimensional  bubble  period  T* 
(where  7*  is  the  first  bubble  period  in  units  of  free-field  bubble  periods)  is  plotted  against  the 
inverse  of  the  non-dimensional  standoff  distance  h*.  For  convenience,  the  results  from 
analyses  with  rigid  wall  boundaries  are  shown  on  the  same  graph  with  results  from  analyses 
with  constant  pressure  boundaries.  The  inverse  standoff  distance  with  an  ordinate  value  of 
0.0  represents  the  situation  in  which  the  boundary  is  located  at  infinity  (i.e.  the  free-field 
case).  Points  to  the  right  of  this  represent  non-dimensional  inverse  standoff  distances  to  a 
rigid  wall,  and  points  to  the  left  are  non-dimensional  inverse  standoff  distances  to  a  constant 
pressure  boundary. 

Also  plotted  in  Figure  4-12,  for  comparison  purposes,  is  the  non-dimensional  period 
versus  the  inverse  of  the  non-dimensional  standoff  distance  (T*  vsersus  (1 //**))  relationship 
predicted  by  an  approximate  analysis  due  to  Herring  (1950).  In  Herring's  analysis,  the  fluid 
is  assumed  to  be  incompressible,  the  gas  bubble  is  assumed  to  have  negligible  internal  energy 
and  remain  nearly  spherical,  and  the  bubble  is  assumed  to  be  fairly  remote  from  the  boundary 
(terms  of  order  l//?*2and  higher  are  neglected).  Herring's  formula  for  the  bubble  period  can 
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Figure  4-12.    Bubble  Period  Variation  with  Standoff  Distance  to  Rigid  and  Constant 
Pressure  Boundaries 
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be  expressed  as 


T*= 


(      (R    IR     )' 

(4-1) 


,   ,    ■     ave      max' 


I 


where  R^e  is  the  average  radius  of  the  free-field  bubble  during  the  first  oscillation,  and  Rmax 
is  the  maximum  free-field  radius.  In  this  equation,  the  upper  (plus)  sign  is  for  a  rigid  wall 
boundary,  and  the  lower  (minus)  sign  is  for  a  constant  pressure  boundary.  The  R^  value  is 
determined  by  integrating  the  free-field  radius  vs  time  curve  over  the  first  period,  and  then 
dividing  by  that  period;  here,  the  free-field  radius  vs  time  curve  from  the  numerical  analysis 
was  used,  and  integrated  numerically. 

It  is  interesting  that  the  T*  versus  \lh*  curve  generated  from  the  numerical  analyses 
has  about  the  same  slope  as  that  predicted  by  Herring's  analysis,  for  the  rigid  boundary  type. 
Herring's  analysis  neglected  terms  of  order  (l//?*2)  and  higher,  and  assumed  small  bubble 
shape  deviations  from  spherical,  so  it  is  apparent  that  equation  (4-1)  could  be  in  error  for 
charges  fairly  near  a  boundary  surface.  Since  Herring  neglected  compressibility  of  the 
surrounding  fluid,  and  changes  in  the  surrounding  flow  field  caused  by  passage  of  the  initial 
shock  wave  (which,  as  seen  in  the  results  in  the  last  chapter,  can  extend  over  an  appreciable 
distance  from  the  bubble),  equation  (4-1)  can  also  be  shown  to  be  in  error  for  charges  very 
remote  from  a  boundary.  In  a  compressible  fluid,  there  is  a  finite  standoff  distance  beyond 
which  the  absence  or  presence  of  a  boundary  can  have  no  affect  upon  the  first  oscillation 
period  of  the  bubble,  owing  to  the  finite  wave  speed  in  a  compressible  medium.  Thus,  in 
Figure  4-12,  for  a  compressible  fluid  there  must  exist  a  finite  region  of  abscissa  values  on 
either  side  of  the  0.0  \lh*  value  for  which  T*  must  be  equal  to  1.0.  This  distance  can  be 
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estimated,  assuming  a  constant  wave  speed  of  about  1500  m/s  in  seawater,  as  half  the  distance 
traveled  by  Shockwave  in  14.47  msec  (the  free-field  bubble  period),  i.e.  about  10.85  meters. 
This  corresponds  to  an  h*  value  of  about  15.4  (and  a  l/h*  value  of  about  0.065). 
Furthermore,  if  the  boundary  is  located  somewhat  closer  it  should  have  still  have  little  effect 
upon  the  bubble  period,  as  most  of  the  oscillation  of  the  bubble  will  have  already  occurred 
before  the  reflection  of  the  primary  shock  wave  gets  back  to  the  bubble. 

Thus,  although  equation  (4-1)  is  the  classical  correction  for  the  effect  of  rigid 
boundaries  and  free-surfaces  on  the  oscillation  period  of  an  explosion  gas  bubble,  a  better 
correction  can  be  made  by  including  fluid  compressibility.  One  question  that  arises  is  how  the 
results  from  this  investigation  this  might  be  applied  underwater  explosion  gas  bubbles 
resulting  from  other  charge  types,  weights,  and  depths.  One  solution  is  to  use  a  simplified  fit 
through  the  data  points  in  Figure  4-12.  One  proposed  fit  is  a  "bilinear"  one  through  the  data 
points  for  either  type  of  boundary.  The  term  "bilinear"  is  meant  to  imply  that,  for  either  a 
rigid  wall  or  a  constant  pressure  boundary  individually,  the  curve  consists  of  two  linear 
segments  (one  of  which  has  zero  slope).  This  proposed  fit  is  shown  as  a  solid  line  in  Figure 
4-12,  and  labeled  as  a  "Bilinear  Approximation."  One  limitation  of  this  fit  is  that  its  accuracy 
decreases  rapidly  for  standoff  distances  of  less  than  two  maximum  free-field  radii,  for  the 
constant  pressure  type  of  boundary  (i.e.  for  \/h*  values  greater  than  0.5).  For  the  rigid  wall 
type  of  boundary,  it  appears  to  be  reasonably  accurate  over  the  entire  range  of  standoff 
distances  shown  in  Figure  4-12. 

The  advantage  of  this  bilinear  approximation  lies  in  the  relative  ease  with  which  it  can 
be  utilized  for  other  charge  type,  weight,  and  depth  parameters,  as  will  be  discussed  shortly. 
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The  functional  form  for  this  proposed  period  correction  factor  for  the  effect  of  a  boundary 
surface  on  the  first  period  of  an  explosion  gas  bubble  is  given  by  the  expression 


T*=l±u 


1    i  \R-JR 


\Vj 


avt      max 


(4-2) 


where  u  is  a  unit  step  function  (zero  if  its  argument  is  negative,  one  if  its  argument  is 
positive),  R^  is  the  average  free-field  bubble  radius  over  the  first  period,  Rmax  is  the  maximum 
free-field  bubble  radius,  and  hT*  is  defined  by 


1       1       1 

TTV^  (4-3) 

ri-      n      w. 


where  h0*  is  an  effective  maximum  standoff  distance  (in  units  of  maximum  free-field  radii)  for 
which  the  period  of  the  bubble  starts  to  be  effected  appreciably  by  the  presence  of  a  boundary 
surface.  In  Figure  4-12,  this  is  the  point  at  which  the  slope  of  the  bilinear  approximation  curve 
changes  from  0.00  to  0.20.  For  the  bilinear  approximation  curve  shown  in  Figure  4-12,  this 
effective  maximum  standoff  distance  is  given  by  \lh0*  =  0.27,  i.e.  h0*  has  a  value  of  about  3.7 
maximum  free-field  radii. 

In  order  to  apply  Herring's  correction  factor,  equation  (4-1),  it  was  necessary  to  know 
the  quantity  R^.  This  implies  knowledge  of  the  free-field  radius  time  history  curve,  but  since 
the  shape  of  this  curve  varies  little  with  the  particular  explosion  parameters,  the  value  of  R^ 
will  generally  be  close  to  0.8  Rmax.  To  apply  equation  (4-2),  it  is  necessary  to  know,  in 
addition,  the  value  of  the  parameter  h0*.  In  general,  this  parameter  may  be  a  function  of  the 
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charge  type,  weight,  and  depth.  However,  it  should  be  relatively  easy  to  determine  an 
empirical  expression  for  this  parameter  for  different  types  of  explosives,  since  for  a  particular 
charge  weight  and  depth,  this  parameter  can  be  determined  by  a  single  accurate  measurement 
of  how  much  Herring's  correction  factor  is  in  error  for  a  single  (well  chosen)  standoff 
distance. 

An  example  should  clarify  this.  Suppose  that  an  experiment  were  conducted  using  a 
particular  charge  type,  weight,  and  depth,  at  a  standoff  distance  of  h*=2.00  from  a  rigid 
boundary,  and  Herring's  correction  factor  was  determined  to  predict  a  period  increase  of 
10%.  The  value  of  RaJRmax  is  then  0.8,  from  equation  (4-1).  If  the  actual  period  increase  for 
this  experiment  was  measured  to  be  4.6%,  then,  from  equation  (4-2),  \/hT*=0.23.  Thus,  from 
equation  (4-3),  \/h0*  =  0.27. 

Again,  the  range  of  validity  of  equation  (4-2)  should  be  kept  in  mind.  From  Figure 
4-12,  an  appropriate  limitation  would  appear  to  be  h*>2.0  for  a  constant  pressure  boundary, 
and  h*>\  .0  for  a  rigid  boundary. 

For  calculation  of  the  first  bubble  period  at  standoff  distances  less  than  two  maximum 
free  field  bubble  radii  from  a  constant  pressure  boundary,  a  better  approximation  is  needed. 
A  "Shifted  Quadratic"  fit  for  the  constant  pressure  boundary  data  shown  in  Figure  4-12  is 
defined  by 


/        \ 


T*=l-u 


\  ^TCPJ 


mTCP 


[hjCP) 


(4-4) 


where  u  is  a  unit  step  function  and 
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1      1      1 

(4-5) 


"TCP  "OCP 


This  characterization  depends  upon  two  parameters.  The  fit  shown  did  not  require  a  linear 
term  in  the  shifted  non-dimensional  inverse  standoff  distance.  The  values  for  the  shifted 
quadratic  fit  shown  in  Figure  4-12  fitting  the  numerical  analysis  data  were  mTCP=0A  and 
tfocjFO.ll.  The  range  of  validity  for  the  fit  defined  by  equations  (4-4)  and  (4-5)  is  A*>1,  i.e. 
it  is  valid  for  any  standoff  distance  from  a  constant  pressure  boundary.  The  fact  that  the  non- 
dimensional  period  change  seems  to  vary  linearly  with  the  non-dimensional  inverse  standoff 
distance  for  a  rigid  wall  boundary  and  as  the  square  of  the  non-dimensional  inverse  standoff 
distance  for  a  constant  pressure  boundary  may  well  have  some  physical  cause  or  explanation. 
2.  The  Effect  of  Boundaries  on  Bubble  Migration  in  a  Compressible  Fluid 

The  bubble  center  of  mass  displacement  time  history  for  the  analyses  in  which  a  rigid 
boundary  and  a  constant  pressure  boundary  were  located  at  an  initial  non-dimensional 
standoff  distance  of  h*=1.380  is  shown  in  Figure  4-13.  This  is  the  same  case  for  which  the 
bubble  volume  equivalent  radius  time  history  curves  were  shown  in  Figure  4-11.  In  Figure 
4-13,  the  displacement  is  taken  to  be  positive  if  it  is  towards  the  boundary,  and  negative  if  it 
is  away  from  the  boundary.  For  the  rigid  wall  boundary  analysis,  this  figure  shows  a  small 
initial  migration  of  the  bubble  away  from  the  boundary  at  early  times  (when  the  bubble  is 
expanding).  At  later  times  (when  the  bubble  is  contracting)  there  is  a  significant  migration 
towards  the  boundary.  This  migration  is  much  larger  than  the  initial  migration  away  from  the 
boundary,  and  the  rate  of  migration  increases  as  the  bubble  contracts  to  its  minimum  volume. 
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Figure  4-13.  Bubble  Center  of  Mass  Displacement  Time  History  for  h*=l  .380  Analyses 
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The  opposite  affect  is  seen  when  the  boundary  is  a  constant  pressure  boundary;  the 
bubble  initially  moves  towards  the  boundary,  and  later  moves  rapidly  away  from  the 
boundary.  The  influence  of  fluid  compressibility  is  also  seen  in  Figures  4-1 1  and  4-13.  There 
is  an  initial  interval  of  time  during  which  neither  the  type  nor  even  the  presence  of  the 
boundary  affects  the  volume  of  the  bubble,  and  it  is  not  migrating. 

The  displacement  time  histories,  whose  overall  characteristic  is  described  above,  also 
show  a  small  amplitude  higher  frequency  component.  The  period  of  this  small  amplitude 
component  is  about  twice  the  time  interval  between  the  start  of  the  analyses  and  the  time 
when  the  bubble  begins  migrating,  which  seems  to  be  roughly  proportional  to  the  standoff 
distance.  Hence,  this  higher  frequency  component  would  appear  to  result  from  the  reflection 
of  the  primary  shock  wave  from  the  boundary. 

The  general  shape  of  the  bubble  displacement  curves,  as  seen  in  Figure  4-13,  has  been 

qualitatively  explained  by  Cole  (1948,  pp.  331-332): 

In  the  case  of  a  rigid  surface,  the  presence  of  the  boundary  interferes  with 
radial  fluid  flow  of  water,  whether  outward  or  inward,  near  a  spherical  surface  in 
its  vicinity.  Initially,  when  the  pressure  in  the  gas  is  in  excess  of  the  hydrostatic 
pressure,  the  water  on  the  side  of  the  bubble  surface  near  the  wall  is  less  readily 
displaced,  and  the  bubble  surface  moves  away  from  the  wall.  The  effect  is 
relatively  small,  however,  because  the  net  pressure  (in  excess  of  hydrostatic)  is 
positive  for  a  short  part  of  the  bubble  period,  and  the  bubble  is  small  during  this 
time.  When  the  pressure  falls  below  hydrostatic,  acceleration  of  flow  toward  the 
bubble  surface  does  not  occur  as  readily  on  the  side  toward  the  wall,  and  the  flow 
must  be  such  as  to  bring  the  surface  nearer  to  the  wall.  A  considerable  amount 
of  momentum  is  imparted  to  a  large  mass  of  water  in  this  way  when  the  bubble  is 
large.  As  the  bubble  contracts,  the  momentum  acquired  becomes  concentrated 
in  a  smaller  mass  of  water  near  the  bubble,  and  the  velocity  of  flow  in  this  region 
increases.  The  bubble  surface  must  then  move  toward  the  wall  with  increasing 
speed  as  if  attracted  to  it.  This  effect  is  so  much  larger  than  the  repulsion  when 
the  pressure  exceeds  hydrostatic  that  the  dominant  motion  is  an  apparent 
attraction  increasing  the  bubble  velocity  toward  the  wall  as  it  contracts,  even 
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though  the  momentum  of  the  flow  is  decreasing  in  the  most  contracted  stages. 

A  free  surface  has  the  opposite  effect  on  bubble  migration,  as  in  this  case 
the  water  at  the  surface  is  free  to  move  but  must  due  so  in  such  a  way  as  to 
equalize  the  pressure  with  that  of  the  atmosphere. 


Herring  (1950)  conducted  an  approximate  analysis  for  the  migration  of  bubbles  near  plane 
rigid  and  free  surfaces,  in  which  the  method  of  images  is  employed  for  motion  of  a  bubble  in  an 
incompressible  fluid,  the  internal  energy  of  the  bubble  is  neglected,  the  bubble  is  assumed  to  be 
not  too  near  the  boundary  surface  (terms  of  order  l/h*2  and  higher  are  neglected)  and  remain 
nearly  spherical,  and  migration  of  the  bubble  is  treated  as  a  small  correction  to  the  motion  of  the 
bubble.  Herring's  formula  for  the  migration  velocity  of  the  center  of  the  bubble  is 


JL*£*-^f'r4*)2«  (4-6) 


where  t  is  the  time,  h  is  the  standoff  distance  between  the  boundary  and  the  center  of  the  bubble, 
and  r  is  the  radius  of  the  bubble.  In  this  equation  the  plus  sign  is  for  a  rigid  boundary,  and  the 
minus  sign  is  for  a  free  surface. 

Together  with  a  knowledge  of  the  function  r(t),  equation  (4-6)  can  be  separated  and 
integrated  to  give  an  equation  for  the  displacement  of  the  center  of  the  bubble.  However,  results 
obtained  from  doing  so  are  known  to  be  in  poor  agreement  with  experimental  results,  significantly 
over  predicting  the  migration  of  the  bubble  late  in  the  collapse  phase  (Cole,  1948,  p.  348).  This 
is  when  the  assumptions  used  by  Herring  in  his  analysis  are  least  accurate,  as  the  migration  of  the 
bubble  at  this  time  is  not  a  small  perturbation  of  the  free-field  bubble  motion,  and  the  time  rate 
of  change  of  the  bubble  radius  is  over  predicted  if  internal  energy  is  neglected. 

Nevertheless,  it  is  worth  noting  that  the  first  term  on  the  right  hand  side  of  equation  (4-6) 
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is  a  periodic  term  which  changes  sign  when  drldt  changes  sign,  i.e.  when  the  bubble  is  at  its 
maximum  radius.  The  second  term  is  monotonic,  and  builds  up  to  a  large  value  when  the 
magnitude  of  drldt  becomes  large.  This  is  in  qualitative  agreement  with  the  behavior  seen  in  the 
displacement  versus  time  results  for  the  various  numerical  analyses. 

As  discussed  earlier,  the  location  of  an  explosion  gas  bubble  at  the  end  of  its  first 
oscillation  is  important  because  this  is  where  the  bubble  pulse  originates  from.  Taking  a  cue  from 
Figure  4-12,  values  for  the  bubble  displacement  at  the  time  the  bubble  reached  minimum  volume 
(which  varies  depending  upon  the  boundary  type  and  initial  standoff  distance)  were  compiled  for 
the  different  numerical  analyses,  and  plotted  as  a  function  of  the  inverse  standoff  distance.  This 
plot  is  shown  in  Figure  4-14,  in  which  the  abscissa  values  represent  the  non-dimensional  inverse 
standoff  distance,  with  values  to  the  right  of  0.0  being  for  a  rigid  boundary  and  values  to  the  left 
of  this  point  being  for  a  constant  pressure  surface,  as  in  Figure  4-12.  The  total  displacements  at 
the  time  of  the  first  minimum  is  plotted  in  units  of  maximum  free-field  radii,  so  that  a  value  of  0.4, 
for  example,  represents  a  center  of  mass  displacement  of  0.4*70.65  cm  from  the  initial  center  of 
the  charge. 

The  effect  of  a  finite  wave  speed  in  a  compressible  fluid  is  again  seen  in  Figure  4-14; 
remote  boundaries  cause  little  or  no  displacement  of  the  bubble.  This  plot  shows  a  general 
linearity  for  bubbles  fairly  near  either  type  of  boundary.  The  only  exception  is  the  case  in  which 
the  bubble  was  initially  closest  to  a  constant  pressure  boundary,  for  which  the  initial  migration 
towards  the  boundary  and  the  increase  in  the  bubble  radius  over  the  free-field  value  cause  the 
bubble  to  actually  contact  the  boundary  when  the  bubble  was  near  maximum  radius.  Neglecting 
this  left  most  point  in  Figure  4-14,  the  results  show  that  the  first  period  displacement  near  either 
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Figure  4-14.    Bubble  CM.  Displacement  at  First  Minimum  Variation  with  Standoff  Distance 
to  Rigid  and  Constant  Pressure  Boundaries 
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type  of  surface  can  be  well  approximated  by  a  bilinear  curve  when  plotted  against  the  inverse 
standoff  distance.  This  curve  is  shown  as  a  dashed  line  in  Figure  4-14.  The  slope  of  the  non-zero 
portion  of  this  curve  is  about  the  same  for  either  type  of  boundary,  again  suggesting  a  genera! 
semi-empirical  formula  for  bubble  migration  near  either  type  of  boundary. 
The  proposed  functional  form  for  this  formula  is 


A7*=« 


/ 

1 


\    *  / 


7 


m7 


(4-7) 


'7 


where  AY*  is  the  displacement  at  the  first  minimum  in  free-field  radii,  u  is  a  unit  step  function,  mY 
is  the  slope  of  the  non-zero  slope  portion  of  the  "bilinear  approximation"  curve,  as  shown  in 
Figure  4-14,  and  hr*  is  a  effective  standoff  distance  for  bubble  migration,  given  by 

1      1       1 

where  h*  is  the  actual  standoff  distance  in  maximum  free-field  bubble  radii,  and  h}*  characterizes 
a  maximum  effective  standoff  distance  for  which  appreciable  bubble  migration  occurs. 

For  the  bilinear  approximation  curve  shown  in  Figure  4-14,  mY  and  h,*  have  the  values 
of  about  0.87  and  6.2,  respectively.  Again,  these  parameters  might  vary  with  charge  type,  weight, 
and  depth,  but  can  be  determined  empirically  without  great  difficulty.  The  proposed  bilinear  fit 
given  by  equations  (4-7)  and  (4-8)  appears,  from  Figure  4-14,  to  be  a  reasonably  good 
approximation  for  total  bubble  migration  at  the  time  of  the  first  minimum  for  standoff  distance 
values  h*>\25  for  the  constant  pressure  type  boundary,  and  h*>\  for  the  rigid  wall  type 
boundary. 
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3.         Bubble  Shape  Departure  From  Spherical 

While  the  time  at  which  the  bubble  reaches  its  minimum  volume  and  the  displacement  of 
the  bubble  at  this  time,  and  the  effect  of  standoff  distance  on  these  quantities,  were  the  primary 
items  of  interest  in  these  analyses,  the  overall  shape  variation  of  an  explosion  gas  bubble 
generated  near  a  boundary  surface  is  also  of  interest.  In  approximate  analytical  derivations  such 
as  for  equation  (4-1),  the  bubble  is  assumed  to  remain  nearly  spherical.  It  is  of  some  interest  to 
see  how  nearly  true  this  assumption  is,  if  only  to  form  some  opinion  as  to  the  validity  of  the 
derivation.  While  the  analysis  code  used  in  these  numerical  analyses  does  not  actually  keep  track 
of  the  location  of  the  boundary  surface  between  different  materials  within  the  mesh,  it  is  possible 
to  get  an  overall  view  of  the  shape  of  the  bubble  by  examining  the  density  of  the  various  mesh 
volumes. 

Figure  4- 1 5  shows  the  results  for  the  analysis  case  in  which  a  plane  constant  pressure 
boundary  was  initially  located  at  a  standoff  distance  of  1.380  maximum  free-field  radii  from  the 
center  of  the  charge.  In  this  figure  the  dark  areas  represent  the  bubble,  and  the  times  are  given 
in  terms  of  the  first  period  of  the  bubble.  The  boundary  plane  is  that  plane  below  which  times  are 
indicated.  The  bubble  is  seen  to  be  fairly  spherical  when  near  its  maximum  volume,  but  to  have 
a  pronounced  "kidney  shape"  by  time  80%  of  its  first  oscillation  period  has  elapsed.  It  then 
continues  deforming,  eventually  assuming  a  "spherical  cap"  shape  by  the  time  it  reaches  minimum 
volume.  This  figure  also  shows  the  migration  of  the  center  of  mass  of  the  bubble  away  from  the 
boundary. 

Figure  4-16  shows  the  variation  of  the  bubble  shape  with  time  when  a  rigid  boundary  was 
located  at  1.380  maximum  free-field  radii  from  the  center  of  the  charge.  Again,  the  dark  area 
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Figure  4-15.    Bubble  Shape  Variation  with  Time  for  Bubble  Located  at  h*=1.380  from  a 
Constant  Pressure  Boundary 
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Figure  4-16.    Bubble  Shape  Variation  with  Time  for  Bubble  Located  at  h*=1.380  from  a 
Rigid  Boundary 
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represents  the  bubble,  and  the  boundary  plane  is  that  plane  below  which  times  are  indicated.  In 
this  figure,  bubble  shapes  are  only  shown  for  the  last  5%  of  the  first  oscillation  period;  at  95%  of 
the  first  bubble  oscillation  period  the  bubble  is  still  nearly  spherical.  By  98%  of  the  first 
oscillation  period,  the  bubble  has  become  "kidney"-shaped. 

From  these  figures,  it  might  be  expected  that  derivations  assuming  bubble  sphericity 
would  be  better  at  this  standoff  distance  for  the  rigid  boundary  case  than  for  the  constant  pressure 
boundary  case.  Examining  the  bilinear  approximation  in  Figure  4-12,  which  can  be  thought  of 
as  a  compressibility  corrected  approximation  to  Herring's  derivation  of  equation  (4-1),  this  would 
seem  to  be  the  case. 
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V.  CONCLUSION 

The  current  understanding  of  the  effects  of  fluid  compressibility  on  the  dynamic  behavior 
of  underwater  explosion  gas  bubbles  is  very  limited;  the  existing  theories  for  underwater 
explosion  gas  bubble  phenomena  neglect  entirely  the  effects  of  fluid  compressibility.  As  a 
consequence,  these  theories  fail  to  accurately  predict  important  phenomena  associated  underwater 
explosion  gas  bubbles.  It  was  for  this  reason  that  the  research  described  in  this  dissertation  was 
undertaken. 

In  the  first  part  of  this  investigation,  it  was  shown  that  a  finite  volume  based  numerical 
analysis  technique  which  incorporated  fluid  compressibility  could  accurately  account  for  the 
behavior  of  an  explosion  gas  bubble  in  the  free-field  case  up  until  the  time  of  the  first  bubble 
minimum,  and  it  was  verified  that  the  radiation  of  energy  in  the  bubble  pulses  does  not  account 
for  all  of  the  energy  loss  that  occurs  between  bubble  oscillation  cycles. 

In  the  second  part  of  this  investigation,  the  behavior  of  underwater  explosion  gas  bubbles 
near  rigid  and  constant  pressure  boundary  surfaces  was  investigated,  including  effects  arising 
because  of  fluid  compressibility.  The  failure  to  account  for  fluid  compressibility  in  the  current 
theoretical  development  has  made  predictions  for  bubble  behavior  based  upon  this  theory  more 
qualitatively  correct  than  quantitatively  accurate.  It  was  found  possible,  when  fluid 
compressibility  is  taken  into  account,  to  quantitatively  describe  the  behavior  of  explosion  gas 
bubbles  quite  close  to  a  boundary  surface. 

By  conducting  analyses  at  a  number  of  different  standoff  distances  from  these  boundaries, 
it  was  discovered  that  first  bubble  oscillation  period  and  the  migration  of  the  bubble  at  the  time 
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it  reached  minimum  volume  could  be  quantitatively  described  over  all  standoff  distances  of 
interest  when  fluid  compressibility  was  taken  into  account.  Functional  relationships  between 
standoff  distance  and  the  first  bubble  period  and  the  net  bubble  migration  were  developed  which 
have  a  form  very  suitable  for  semi-empirically  scaling  the  results  of  this  investigation  to  situations 
involving  other  charge  types,  weights,  and  depths. 

It  is  recommended  that  a  series  of  carefully  designed  experiments  be  conducted  to  both 
verify  the  formulas  developed  in  this  investigation  and  to  develop  the  neccessary  empirical 
constants  used  in  these  formulas,  so  that  the  effects  of  boundary  surfaces  can  be  better  taken  into 
account  in  the  future.  These  results  of  this  investigation  should  form  a  foundation  for  improving 
the  accuracy  of  calculations  for  the  whipping  response  of  ships  and  shipboard  equipment  due  to 
detonation  of  nearby  mines  or  torpedoes. 
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